# Inferring synergistic and antagonistic interactions in mixtures of exposures

Shounak Chattopadhyay<EMAIL_ADDRESS>Department of Statistical
Science, Duke University  Stephanie M. Engel<EMAIL_ADDRESS>Department of Epidemiology, Gillings School of Global Public Health,
University of North Carolina at Chapel Hill  David Dunson<EMAIL_ADDRESS>Department of Statistical Science, Duke University

###### Abstract

There is abundant interest in assessing the joint effects of multiple
exposures on human health. This is often referred to as the mixtures problem
in environmental epidemiology and toxicology. Classically, studies have
examined the adverse health effects of different chemicals one at a time, but
there is concern that certain chemicals may act together to amplify each
other’s effects. Such amplification is referred to as synergistic interaction,
while chemicals that inhibit each other’s effects have antagonistic
interactions. Current approaches for assessing the health effects of chemical
mixtures do not explicitly consider synergy or antagonism in the modeling,
instead focusing on either parametric or unconstrained nonparametric dose
response surface modeling. The parametric case can be too inflexible, while
nonparametric methods face a curse of dimensionality that leads to overly
wiggly and uninterpretable surface estimates. We propose a Bayesian approach
that decomposes the response surface into additive main effects and pairwise
interaction effects, and then detects synergistic and antagonistic
interactions. Variable selection decisions for each interaction component are
also provided. This Synergistic Antagonistic Interaction Detection (SAID)
framework is evaluated relative to existing approaches using simulation
experiments and an application to data from NHANES.

Keywords: Additive model, Bayesian, Dose response surface, Nonparametric
Regression, Pairwise interaction, Positivity constraint, synergy, Exposure
mixture

##  1 Introduction

There is considerable concern that humans are exposed to a variety of
potentially adverse chemicals, and these exposures may have adverse health
effects. Classically, such health effects have been assessed through one
exposure at a time studies, either collecting in vitro or in vivo data at
different doses of a single chemical or focusing analyses of observational
epidemiology studies on single exposures. Then, in order to predict the
overall health effect of a mixture of different exposures, one needs to make
strong assumptions, such as additivity. Unfortunately, such predictions will
misestimate an individual’s true adverse health risk if certain chemicals
interact. Of particular concern are synergistic interactions in which the
adverse effect of one chemical is increased due to the presence of another
chemical. If regulatory agencies are unaware of such synergistic interactions
in setting pollution guidelines, they may inadvertently permit substantial
pollution-induced mortality and morbidity. On the other hand, antagonistic
interactions may lead to additive models over-predicting risk in which case
certain chemicals may be over-regulated.

In this paper, our goal is to analyze the effect of heavy metal exposure on
human kidney function. The impact of metal exposures on human renal function
has garnered considerable attention from the epidemiological community. A
review of the existing literature shows evidence of degraded kidney function
following prolonged exposure to heavy metals (Pollack et al.,, 2015; Luo and
Hendryx,, 2020). Most of this literature is based on simple one exposure-at-a-
time correlation analyses from observational epidemiology data. It is clearly
of interest to study joint effects of multiple metals, while adjusting for
covariates that may act as potential confounding variables. Mechanistically,
it seems unlikely that metals have a simple additive effect on kidney
function. For example, healthy kidney function may continue with a single
metal exposure at relatively low doses, but as dose increases and/or
additional metal exposures are added it is likely that kidney function may
worsen rapidly. Such a dose response surface would be reflective of a
synergistic interaction. In our analyses, we focus on 2015-16 data from
NHANES. These data contain information on heavy metals found in spot urine
collections and also on urine creatinine, which can be used as a marker of
kidney function when urinary dilution issues are appropriately accounted for.

There is a clear need for new statistical tools for accommodating interactions
in assessing the health effects of mixtures of chemical exposures. For a
review of recent developments, refer to Joubert et al., (2022). A broad set of
strategies have been taken in this literature. The first extends linear
regression to include a quadratic term characterizing pairwise interactions;
for recent examples, refer to Ferrari and Dunson, (2021); Wang et al., (2019).
This approach can trivially identify synergistic versus antagonistic
interactions through the signs of the quadratic coefficients, but relies on a
highly restrictive parametric model. An alternative is to rely on
nonparametric regression; for example, using Bayesian Additive Regression
Trees (BART) (Chipman et al.,, 2010) or Gaussian processes (GP) (Williams and
Rasmussen,, 2006). Bayesian Kernel Machine Regression (BKMR) (Bobb et al.,,
2015) implements GP regression with variable selection motivated by the
mixtures problem. These approaches do not identify synergistic or antagonistic
interactions, and can lead to overly-wiggly and difficult to interpret dose
response surfaces. MixSelect (Ferrari and Dunson,, 2020) bridges between
parametric and nonparametric approaches via an additive expansion. However, it
lacks flexibility to characterize nonlinear main effects and interactions due
to its reliance on quadratic regression. To simplify dose response modeling,
Molitor et al., (2010) propose profile regression based on clustering
exposures, while Czarnota et al., (2015) use a single index model based on a
weighted sum of the exposures. Neither approach allows for inferences on
interactions.

We propose a nonparametric Bayesian approach for identifying synergistic or
antagonistic interactions between $p$ exposures. The dose response surface is
decomposed additively into $p$ main effects and $\binom{p}{2}$ pairwise
interactions. Ruling out higher order interactions substantially reduces
dimensionality, while aiding interpretability. Each pairwise interaction is
decomposed into the difference of two non-negative functions, facilitating
selection of synergistic or antagonistic effects. The proposed synergistic
antagonistic interaction detection (SAID) approach is shown to improve over
unconstrained nonparametric regression, leading to comparatively higher
statistical power in simulation studies. To avoid modeling complicated
pairwise surfaces, we focus on a special case to dramatically reduce
dimensionality and obtain substantial computational gains. The resulting
interaction surface retains a high degree of flexibility in capturing
nonlinear surfaces. We also outline an approach to carry out variable
selection on the pairwise interactions. The SAID approach employs a
computationally efficient Markov chain Monte Carlo (MCMC) algorithm to sample
from the posterior distribution of the unknown quantities.

The paper is outlined as follows. We describe the urine creatinine and metal
exposure data in Section 2 and outline challenges involved with analyzing
these data. Section 3 provides details of our SAID approach. Section 4
compares SAID and existing competitors in simulation experiments, highlighting
benefits of using SAID in terms of estimation accuracy, valid uncertainty
quantification, and variable selection. Section 5 uses SAID to investigate
synergistic and antagonistic interactions of metal exposures in predicting
urine creatinine levels based on NHANES 2015-16 data. Section 6 includes
concluding remarks.

##  2 Kidney Function Data Analysis

###  2.1 Motivation

We are interested in assessing the impact of exposure to heavy metals on human
renal function. The concentration of creatinine in the blood or urine is an
established biomarker of kidney function (Kashani et al.,, 2020; Barr et al.,,
2005). Exposure to heavy metals has been linked to changes in renal function
and kidney damage (Pollack et al.,, 2015). In addition to being indicative of
renal function, creatinine levels are also related to muscle mass (Forbes and
Bruining,, 1976; Baxmann et al.,, 2008). Urine creatinine has been shown to be
positively associated with serum creatinine levels in studies excluding
individuals with chronic kidney disease (CKD) (Jain,, 2016). Existing studies,
such as the one in Luo and Hendryx, (2020), show a statistical association
between biomarkers of chronic kidney disease and blood levels of the heavy
metals cobalt (Co), chromium (Cr), mercury (Hg), and lead (Pb) using NHANES
data. Kim et al., (2015) find an association between CKD and elevated levels
of cadmium (Cd) in blood. These studies primarily focus on marginal
associations between individual chemicals and CKD. Our focus is instead on
studying how multiple heavy metal exposures relate to kidney function measured
through urine creatinine, with an emphasis on identifying synergistic and
antagonistic interactions.

###  2.2 Data Description

We analyze data that were collected by NHANES for the year $2015$. We consider
urine analyte levels of the $13$ heavy metals: Antimony (Sb), Barium (Ba),
Cadmium (Cd), Cesium (Cs), Cobalt (Co), Lead (Pb), Manganese (Mn), Molybdenum
(Mo), Strontium (Sr), Thallium (Tl), Tin (Sn), Tungsten (W), and Uranium (U)
as exposure variables and the level of urine creatinine (uCr) as the response
variable. The unit of measurement for exposure variables is $\mu$g/mL while
the unit for the response variable is mg/dL. We also adjust for age (in
years), sex (male or female), ethnicity (Non-Hispanic White, Non-Hispanic
Black, Mexican American, Other Hispanic, and Other), and body mass index
(BMI).

We start off by only considering the data for $2300$ individuals for which
values of the response variable (uCr) are not missing. In this data set, we
first remove the individuals having serious kidney disease. For this, we used
the albumin-to-creatinine ratio (ACR) and removed the subjects who had
albuminuria, defined as their ACR satisfying $\mbox{ACR}\geq 30$ mg/g. The
resulting data set has $2008$ individuals. Furthermore, the response variable
for one individual was less than the limit of detection (LOD); this singular
data point was removed for purposes of analysis. For the metal exposures,
there are both missing data points and data points which are lower than their
respective LODs. We remove any individuals with an exposure entry missing.
Finally, we removed the entries of the covariate values which were missing. In
summary, neither the response variable nor the covariate variables considered
in our analysis have any missing data or data less than the LOD; however, the
metal exposure data contain points less than the LOD. After removing the
individuals missing response, covariate, and exposure measurements, the final
sample size is $n=1979$.

We look at various summary measures of the response and exposure variables.
The uCr levels are right skewed, with mean $119.67$ mg/dL, median $104.00$
mg/dL, and semi-interquartile range (SIQR) $52.50$ mg/dL. The heavy metal
exposure variables are mild to moderately correlated, with a heat plot of the
correlation matrix shown in Figure 1. More specifically, $58$ out of the $78$
pairwise correlation coefficients are less than $0.20$. The median correlation
is $0.10$ with an SIQR of $0.08$, while the maximum correlation is $0.51$,
between Barium and Strontium. Furthermore, there is wide variation between the
range of different exposure variables. For example, the maximum recorded
exposure to Cesium in the sample is $110.87$ $\mu$g/mL, while that of Uranium
is only $0.77$ $\mu$g/mL.

Figure 1: Heat plot showing correlations between heavy metals in NHANES
2015-16 data.

###  2.3 Urinary Dilution

In NHANES, the urine data is collected from spot collections. The urine sample
collection procedure is regulated to ensure consistent specimen collection
across differing water loadings of the subjects. Regardless, due to the nature
of data collection, the concentrations of the heavy metals and creatinine
levels can vary substantially across the sampled participants due to different
hydration levels. For this reason, we adjust the raw response (creatinine) and
exposure variables (heavy metals) measured in the urine using the urine flow
rate (UFR, measured in mL/min) of the subjects. Adjusting for dilution using
urine flow rate also helps reduce wide variation in exposure measurements.
Adjusting urinary measures for variation in water loading across individuals
using urine flow rates is prevalent in existing literature (Jeng et al.,,
2021; Middleton et al.,, 2016; Hays et al.,, 2015). Suppose for the $i$-th
subject, the urine creatinine level is $C_{i}$ and the urinary concentrations
of the heavy metals is the vector ${\bf
M}_{i}=(M_{i1},\ldots,M_{ip})^{\text{T}}$. Provided the urine flow rate for
the $i$-th subject is $\tau_{i}>0$, we adjust for urinary dilution by
multiplying both the original response and exposure variables by $\tau_{i}$.
That is, the urine flow adjusted response variable and exposure vector are
$\tau_{i}C_{i}$ and $\tau_{i}{\bf M}_{i}$, respectively, which we refer to as
the dilution-adjusted response and exposure variables, respectively.

###  2.4 Issues with existing approaches

As two different types of state-of-the-art methods, we apply MixSelect and
BKMR on the heavy metals and creatinine data. We (natural) log transform the
response and exposures prior to analysis. Both BKMR and MixSelect conduct
Bayesian variable selection, providing posterior inclusion probabilities
(PIPs) for each exposure. Excluding exposures having PIPs $<0.5$, BKMR
excludes Barium (PIP $\approx$ 0), Lead (PIP $\approx$ 0), and Tin (PIP =
0.13), while MixSelect excludes Cobalt (PIP = 0.29) and Uranium (PIP = 0.11).
With these exposures excluded, BKMR produces a $10$ dimensional dose response
surface in the remaining metals. Although inferences on main effects and
pairwise interactions could rely on examining univariate and bivariate cross-
sections of the $10$-dimensional surface, such results are exploratory and
difficult to interpret. In contrast, MixSelect provides PIPs for both the main
effects and pairwise interactions. However, it is not reassuring that
MixSelect excludes different exposures than BKMR. In addition, all the
interaction PIPs for MixSelect are $<0.11$, indicating that none of the
interactions are selected.

We are motivated by these preliminary results to develop an approach that is
flexible enough to characterize nonlinear dose response surfaces, while
allowing us to formally detect pairwise interactions that are synergistic or
antagonistic. Flexible nonparametric dose response surface methods can
characterize synergistic or antagonistic pairwise interactions but they tend
to be hidden within a flexible multivariate surface, as highlighted for BKMR
above. We describe our proposed Synergistic Antagonistic Interaction Detection
(SAID) approach in Section 3. SAID will be applied to the motivating metals
and creatinine data in Section 5.

##  3 Structured Interaction Modeling Approach

###  3.1 Basic Modeling Structure

Denote the health outcome of interest by $y_{i}$, for individuals
$i=1,\ldots,n$; we assume $y_{i}\in\mathbb{R}$ but the approach can be easily
modified to allow $y_{i}\in\\{0,1\\}$ or ordered categorical or count
responses. We denote the exposures as
$\mathbf{x}_{i}=(x_{i1},\ldots,x_{ip})^{\text{T}}\in[0,1]^{p}$, assuming
without loss of generality they each fall within $[0,1]$, and the covariates
as $\mathbf{z}_{i}\in\mathbb{R}^{q}$. The model we consider is given by

$y_{i}=H(\mathbf{x}_{i})+\eta^{\text{T}}\mathbf{z}_{i}+\epsilon_{i},\quad\epsilon_{i}\sim
N(0,\sigma^{2}),$ (1)

where $H(\mathbf{x})=H(x_{1},\ldots,x_{p})$ is the dose response function of
the exposures, and we follow common practice in adjusting for the covariates
linearly. Without loss of generality, we assume that higher values of the
health outcome $y$ represent worse health in interpreting exposure effects. As
in Wei et al., (2020); Brezger and Lang, (2006), we characterize the dose
response function of the exposures $H({\bf x})$ via an additive expansion into
main effects and pairwise interaction terms:

$H(x_{1},\ldots,x_{p})=\alpha+\sum_{j=1}^{p}f_{j}(x_{j})+\sum_{1\leq u<v\leq
p}h_{uv}(x_{u},x_{v}),$ (2)

where $\alpha$ is an intercept, $f_{j}(x_{j})$ is the main effect of the $j$th
exposure for $j=1,\ldots,p$, and $h_{uv}(x_{u},x_{v})$ for $1\leq u<v\leq p$
is a pairwise interaction. The decomposition in (2) allows interpretation of
different components in a factorization of the dose response surface $H$ into
main effects and pairwise interactions.

The primary innovation of this paper is our approach for inferences on the
pairwise interaction component. Section 3.2 introduces the general structure
of our model for the $h_{uv}$s. Section 3.3 describes an approach for
interaction selection. Section 3.4 describes our model for the main effects.
Section 3.5 contains details on posterior computation.

###  3.2 Modeling pairwise interactions

We propose a model for the pairwise interactions $h_{uv}$s. Assume that
$h_{uv}$ is continuous and admits finite partial derivatives up to second
order. For $\mathbf{a}=(a_{1},\ldots,a_{d})^{\text{T}}\in\mathbb{R}^{d}$, we
let $\mathbf{a}^{2}=(a_{1}^{2},\ldots,a_{d}^{2})^{\text{T}}\in\mathbb{R}^{d}$.
For a non-negative function $f:[0,1]^{2}\to[0,\infty)$, we say $f\equiv 0$ if
$f(x_{1},x_{2})=0$ for all $(x_{1},x_{2})\in[0,1]^{2}$; otherwise, if
$f(x_{1},x_{2})>0$ for at least one $(x_{1},x_{2})\in[0,1]^{2}$, $f\not\equiv
0$.

In conducting inferences on interactions in environmental epidemiology, it is
of primary interest to assess whether exposures work together to magnify their
health effects (synergy), tend to block each other’s effects (antagonism), or
have effectively no interaction (null). Hence, for exposures $u$ and $v$, we
focus on the following classes of interactions.

  1. 1.

Synergistic if $h_{uv}(x_{u},x_{v})\geq 0$ for all $(x_{u},x_{v})\in[0,1]^{2}$
with at least one strict inequality.

  2. 2.

Antagonistic if $h_{uv}(x_{u},x_{v})\leq 0$ for all
$(x_{u},x_{v})\in[0,1]^{2}$ with at least one strict inequality.

  3. 3.

Null if $h_{uv}(x_{u},x_{v})=0$ for all $(x_{u},x_{v})\in[0,1]^{2}$.

Our primary interest is in classifying interactions $h_{uv}$ in terms of
Definitions 1-3. Although we will not rule out other types of interaction
surfaces a priori, we develop inference approaches that are targeted towards
this three class hypothesis testing problem. This is accomplished in a
Bayesian manner with a carefully-structured model for $h_{uv}$ that shrinks
towards the space of synergistic, antagonistic, and null interactions.

If we knew a priori  that $h_{uv}$ was either synergistic or null, we could
let $h_{uv}=P_{uv}$, where $P_{uv}\geq 0$. If $P_{uv}(x_{u},x_{v})>0$ for at
least one $(x_{u},x_{v})\in[0,1]^{2}$ then $h_{uv}$ is synergistic; otherwise,
$P_{uv}\equiv 0$ and $h_{uv}$ is a null interaction. Similarly, if we knew a
priori that $h_{uv}$ was either antagonistic or null, we could let
$h_{uv}=-N_{uv}$, where $N_{uv}\geq 0$. However, the sign of the interaction
$h_{uv}$ is usually unknown and pairwise interactions may be only
approximately synergistic or antagonistic. To allow for such complexities, we
let

$h_{uv}(x_{u},x_{v})=P_{uv}(x_{u},x_{v})-N_{uv}(x_{u},x_{v}),$ (3)

for all $(x_{u},x_{v})\in[0,1]^{2}$, where
$P_{uv},N_{uv}:[0,1]^{2}\to[0,\infty)$ are non-negative functions. If
$P_{uv}\not\equiv 0$ and $N_{uv}\equiv 0$, $h_{uv}$ is synergistic; while if
$P_{uv}\equiv 0$ and $N_{uv}\not\equiv 0$, $h_{uv}$ is antagonistic. If both
$P_{uv}=N_{uv}\equiv 0$, $h_{uv}$ is null. All three cases are contained in
the restriction $\int P_{uv}\int N_{uv}=0$, where $\int f$ is shorthand for
$\int_{[0,1]^{2}}f(x_{1},x_{2})\,dx_{1}dx_{2}$ for a bivariate function $f$.
By penalizing $\int P_{uv}\int N_{uv}$, we can favor $h_{uv}$ close to
synergistic, antagonistic, or null.

We decompose $P_{uv}$ and $N_{uv}$ as products of non-negative univariate
functions:

$\displaystyle h_{uv}(x_{u},x_{v})$
$\displaystyle=P_{uv}(x_{u},x_{v})-N_{uv}(x_{u},x_{v})$ (4)
$\displaystyle=P_{uv,1}(x_{u})P_{uv,2}(x_{v})-N_{uv,1}(x_{u})N_{uv,2}(x_{v}).$

This model is much more flexible than commonly used quadratic regression,
which lets $h_{uv}(x_{u},x_{v})=\gamma_{uv}x_{u}x_{v}$. To model
$P_{uv,1},P_{uv,2}$, $N_{uv,1},N_{uv,2}$ as flexible one-dimensional non-
negative functions, we rely on squaring B-spline expansions (De Boor,, 1978)
as follows:

$\displaystyle
P_{uv,1}(x_{u})=\\{\mathbf{s}_{u}(x_{u})^{\text{T}}\theta_{uv,1}\\}^{2}$
$\displaystyle,\quad
P_{uv,2}(x_{v})=\\{\mathbf{s}_{v}(x_{v})^{\text{T}}\phi_{uv,1}\\}^{2}$ (5)
$\displaystyle
N_{uv,1}(x_{u})=\\{\mathbf{s}_{u}(x_{u})^{\text{T}}\theta_{uv,2}\\}^{2}$
$\displaystyle,\quad
N_{uv,2}(x_{v})=\\{\mathbf{s}_{v}(x_{v})^{\text{T}}\phi_{uv,2}\\}^{2},$

where $\mathbf{s}_{u}(x_{u})=(s_{u1}(x_{u}),\ldots,s_{um}(x_{u}))^{\text{T}}$
denote B-splines for the $u$-th variable, chosen so that $s_{uj}(0)=0$ for
$j=1,\ldots,m$. This is obtained by discarding the intercept spline; refer to
Section A of the Supplementary Material for further details. Constraining
$s_{uj}(0)=0$ for $j=1,\ldots,m$ and $u=1,\ldots,p$ ensures that the
interaction $h_{uv}$ satisfies $h_{uv}(x_{u},0)=0$ for all $x_{u}\in[0,1]$ and
$h_{uv}(0,x_{v})=0$ for all $x_{v}\in[0,1]$, thereby making $h_{uv}$
identifiable. For further details on identifiability of pairwise interactions,
we refer the reader to Sections A and C of the Supplementary Material.

To estimate $h_{uv}$ in a Bayesian manner, we define a prior distribution
$\pi(\mathbf{\Psi}_{uv})$ on the basis coefficients
$\mathbf{\Psi}_{uv}=(\theta_{uv,1}^{\text{T}},\phi_{uv,1}^{\text{T}},\theta_{uv,2}^{\text{T}},\phi_{uv,2}^{\text{T}})^{\text{T}}$.
In order to formulate $\pi(\mathbf{\Psi}_{uv})$, we first define a prior
distribution $\pi_{0}(\mathbf{\Psi}_{uv})$ and then augment it with a penalty
term based on $\mathcal{Q}(P_{uv},N_{uv})=\int P_{uv}\int N_{uv}$. The initial
prior $\pi_{0}(\mathbf{\Psi}_{uv})$ is defined as

$\theta_{uv,1},\phi_{uv,1}\sim
N\left(0,\nu^{2}\tau_{uv,1}^{2}\Sigma_{0}\right),\quad\theta_{uv,2},\phi_{uv,2}\sim
N\left(0,\nu^{2}\tau_{uv,2}^{2}\Sigma_{0}\right),\\\ $ (6)

where $\Sigma_{0}$ is a P-spline covariance matrix as in Lang and Brezger,
(2004); for more details, refer to Section B of the Supplementary Material.
Let $w_{uv}=(\tau_{uv,1},\tau_{uv,2})$ denote the $uv$-th interaction specific
variance parameters. Then, the augmented prior is

$\pi(\mathbf{\Psi}_{uv}\mid\kappa_{uv},w_{uv},\nu)\propto\pi_{0}(\mathbf{\Psi}_{uv}\mid
w_{uv},\nu)\exp\\{-\kappa_{uv}\mathcal{Q}(P_{uv},N_{uv})\\}.$

By letting
$A_{u}=\int\mathbf{s}_{u}(x_{u})\mathbf{s}_{u}(x_{u})^{\text{T}}\,dx_{u}$ for
$u=1,\ldots,p$, $\mathcal{Q}(P_{uv},N_{uv})$ can be expressed as

$\mathcal{Q}(P_{uv},N_{uv})=\widetilde{\mathcal{Q}}(\mathbf{\Psi}_{uv})=\left(\theta_{uv,1}^{\text{T}}A_{u}\theta_{uv,1}\right)\left(\phi_{uv,1}^{\text{T}}A_{v}\phi_{uv,1}\right)\left(\theta_{uv,2}^{\text{T}}A_{u}\theta_{uv,2}\right)\left(\phi_{uv,2}^{\text{T}}A_{v}\phi_{uv,2}\right).$

We can therefore rewrite $\pi(\cdot\mid\kappa_{uv},w_{uv},\nu)$ as

$\pi(\mathbf{\Psi}_{uv}\mid\kappa_{uv},w_{uv},\nu)\propto\pi_{0}(\mathbf{\Psi}_{uv}\mid
w_{uv},\nu)\exp\\{-\kappa_{uv}\widetilde{\mathcal{Q}}(\mathbf{\Psi}_{uv})\\}.$
(7)

The prior is completed with hyperpriors for the variance and penalty
parameters:

$\tau_{uv,1},\tau_{uv,2}\sim C^{+}(0,1),\quad\log(\kappa_{uv})\sim
N(0,1),\quad\nu\sim C^{+}(0,1),$ (8)

where $C^{+}(0,1)$ denotes a half-Cauchy prior.

Prior (7)-(8) is chosen to have a global-local shrinkage form motivated by the
horseshoe prior (Carvalho et al.,, 2009). Small values of the global parameter
$\nu$ favor most of the interactions $h_{uv}\approx 0$, while the heavy-tailed
prior on local parameters $\tau_{uv,1}$ and $\tau_{uv,2}$ allow certain
interactions to have $P_{uv}$ and $N_{uv}$ components arbitrarily far from
zero. The role of the penalty is to favor $h_{uv}$s that are close to
synergistic, antagonistic, or null. However, the prior on $\kappa_{uv}$ allows
the penalty to be greatly relaxed for certain pairs of exposures.

To investigate the effect of the penalty parameter $\kappa_{uv}$, we let
$p=2$, corresponding to a single interaction $h_{12}$, and simulate draws from
the prior $\pi$ in (7). We fix $\tau_{12,1}=\tau_{12,2}=\nu=1$ and vary
$\kappa_{12}\in\\{0,0.1,1,10,100,1000\\}$. Given a value of $\kappa_{12}$, we
draw $\mathbf{\Psi}_{12}\sim\pi$ and compute $\mathcal{W}=\\{\int
h_{12}^{+}(x_{u},x_{v})\,dx_{u}dx_{v}\\}\\{\int
h_{12}^{-}(x_{u},x_{v})\,dx_{u}dx_{v}\\}$, where
$h_{12}^{+}(x_{u},x_{v})=\max\\{h_{12}(x_{u},x_{v}),0\\}$ and
$h_{12}^{-}(x_{u},x_{v})=\max\\{-h_{12}(x_{u},x_{v}),0\\}$. It is clear that
$\mathcal{W}=0$ if and only if $h_{12}$ is either synergistic, antagonistic,
or null. In Figure 2, we plot the proportion of prior draws of $\mathcal{W}$
less than $0.001$ for each value of $\kappa_{12}$. The proportion of prior
samples of $\mathcal{W}$ less than $0.001$ increases from $0.13$ to $0.87$ as
we increase the penalty $\kappa_{12}$ from $0$ to $1000$. We also looked at
the maximum of $\mathcal{W}$ over prior draws, which decreased from $939.03$
to $0.05$ as $\kappa_{12}$ increased from $0$ to $1000$. This shows that the
proposed prior distribution is flexible enough to capture arbitrary
interactions for small values of $\kappa_{12}$, while favoring interactions
that are very close to being synergistic, antagonistic, or null for large
values of $\kappa_{12}$.

Figure 2: Plot showing proportion of prior draws of $\mathcal{W}$ less than
$0.001$ as $\kappa_{12}$ increases. $\mathcal{W}$ is a measure of deviation of
the interaction $h_{12}$ from being synergistic, antagonistic, or null.

We also tried other approaches to model $h_{uv}$. Instead of squaring
unconstrained functions, one could potentially use a Bayesian monotone spline
formulation based on restricting the sign of the basis coefficients. However,
we found squaring an unconstrained function to provide better estimates,
particularly when the non-negative function being modeled is near $0$. We also
tried using a tensor product spline model for the bivariate functions, but
using products of univariate functions had better power. We describe these
alternative approaches in more detail in Section D of the Supplementary
Material.

###  3.3 Variable Selection

To select the non-null interactions $h_{uv}$s, we first decompose

$h_{uv}(x_{u},x_{v})=h_{uv}^{+}(x_{u},x_{v})-h_{uv}^{-}(x_{u},x_{v}),$ (9)

where $h_{uv}^{+}(x_{u},x_{v})=\max\\{h_{uv}(x_{u},x_{v}),0\\}\geq 0$ and
$h_{uv}^{-}(x_{u},x_{v})=\max\\{-h_{uv}(x_{u},x_{v}),0\\}\geq 0,$ for any
$(x_{u},x_{v})\in[0,1]^{2}$. The interaction $h_{uv}$ is synergistic if and
only if $h_{uv}^{+}\not\equiv 0$ and $h_{uv}^{-}\equiv 0$, antagonistic if and
only if $h_{uv}^{+}\equiv 0$ and $h_{uv}^{-}\not\equiv 0$, and null if and
only if $h_{uv}^{+}=h_{uv}^{-}\equiv 0$. These conditions may be rewritten in
terms of their integrals since for a continuous non-negative function $P$,
$P\equiv 0$ if and only if $\int P=0$ and $P\not\equiv 0$ if and only if $\int
P>0$.

Because we are using a continuous shrinkage prior that places zero probability
on $h_{uv}^{+}$ or $h_{uv}^{-}$ being exactly zero, we treat the integrals
$\int h_{uv}^{+}$ and $\int h_{uv}^{-}$ as effectively zero if they are below
a cutoff $c_{0}>0$. We call $c_{0}$ the integral cutoff. Let
$C_{uv}=\left\\{\int h_{uv}^{+}\leq c_{0}\right\\}$ and $D_{uv}=\left\\{\int
h_{uv}^{-}\leq c_{0}\right\\}$ denote the events that the sizes of
$h_{uv}^{+}$ and $h_{uv}^{-}$ are less than the integral cutoff $c_{0}$,
respectively. In practice, we compute Monte Carlo estimates of the posterior
probabilities $\mathbf{P}\left(C_{uv}\cap
D_{uv}\right),\mathbf{P}\left(C_{uv}^{c}\cap D_{uv}\right),$ and
$\mathbf{P}\left(C_{uv}\cap D_{uv}^{c}\right)$; these are then used as the
estimates of the posterior probabilities that $h_{uv}$ is null, synergistic,
and antagonistic, respectively. In particular, the posterior inclusion
probability (PIP) of the interaction $h_{uv}$ is given by

$\mbox{PIP}(h_{uv})=1-\mathbf{P}(C_{uv}\cap D_{uv})=\mathbf{P}(C_{uv}^{c}\cup
D_{uv}^{c})$ (10)

In practice, we tried with different integral cutoff values
$c_{0}=0.005,0.01,0.05,0.10$. We provide further details on sensitivity
analysis to $c_{0}$ in Section F of the Supplementary Material. Based on the
sensitivity analysis, we recommend using an integral cutoff of $c_{0}=0.01$,
after standardizing $y$ to have variance $1$ prior to fitting SAID. As a rough
rule of thumb, loosely based on Kass and Raftery, (1995), one can view
interactions having $\mbox{PIP}(h_{uv})\in(0.5,0.75)$ as barely worth a
mention in terms of weight of evidence against the null,
$\mbox{PIP}(h_{uv})\in[0.75,0.95)$ as weakly to moderately suggestive,
$\mbox{PIP}(h_{uv})\in[0.95,0.99)$ as strong evidence, and
$\mbox{PIP}(h_{uv})\geq 0.99$ as very strong evidence.

###  3.4 Main Effects and Other Parameters

For reasons of identifiability, we assume the main effects
$f_{1},\ldots,f_{p}$ either satisfy $\int_{0}^{1}f_{j}(x_{j})\,dx_{j}=0$ for
each $j=1,\ldots,p$ or $f_{j}(0)=0$ for $j=1,\ldots,p$. We call the former set
of conditions as integral constraints and the latter as origin constraints.
Suppose that for exposure $j$, $\\{b_{j,1}(\cdot),\ldots,b_{j,d}(\cdot)\\}$
represent one-dimensional basis functions, such as B-splines. To enforce the
identifiability conditions, we let the functions $b_{j,1},\ldots,b_{j,d}$
either satisfy $\int_{0}^{1}b_{j,u}(x_{j})\,dx_{j}=0$ for $u=1,\ldots,d$,
$j=1,\ldots,p$ if the main effects satisfy the integral constraint, or
$b_{j,u}(0)=0$ for $u=1,\ldots,d$, $j=1,\ldots,p$ if the main effects satisfy
the origin constraint. We provide further details in Sections A and C of the
Supplementary Material. Let
$\mathbf{b}_{j}(x_{j})=(b_{j,1}(x_{j}),\ldots,b_{j,d}(x_{j}))^{\text{T}}$. We
now model $f_{j}$ as

$f_{j}(x_{j})=\sum_{u=1}^{d}b_{j,u}(x_{j})\gamma_{j,u}=\mathbf{b}_{j}(x_{j})^{\text{T}}\,\mathbf{\Gamma}_{j},$
(11)

where $\mathbf{\Gamma}_{j}=(\gamma_{j,1},\ldots,\gamma_{j,d})^{\text{T}}$. To
estimate the coefficient vector $\mathbf{\Gamma}_{j}$, we assume a univariate
P-spline prior distribution (Lang and Brezger,, 2004) on
$\mathbf{\Gamma}_{j}$. With an appropriate positive-definite matrix
$\Sigma_{M}$, we can write this prior as

$\mathbf{\Gamma}_{j}\mid\lambda_{j}\sim
N\left(0,\frac{\Sigma_{M}}{\lambda_{j}}\right),\quad\lambda_{j}\sim G(a,a),$
(12)

where $\lambda_{j}>0$ is a scale parameter corresponding to the $j$-th main
effect. The sensitivity to the choice of hyperparameter $a$ for the
distribution of $\lambda_{j}$ has been investigated in Lang and Brezger,
(2004); we found $a=0.5$ to work the best across simulations and applications.

To complete prior specification for the parameters in (1)-(2), we put vague
prior distributions on the intercept $\alpha$ and covariate effects $\eta$,
and a non-informative prior on the measurement error variance $\sigma^{2}$:

$\displaystyle\alpha$ $\displaystyle\sim N(0,10^{4}),\quad\mathbf{\eta}\sim
N(0,10^{4}\,\mathbb{I}_{q}),\quad\sigma^{2}\sim\sigma^{-2}.$ (13)

###  3.5 Posterior Sampling

We rely on a Hamiltonian Monte Carlo (HMC) (Neal et al.,, 2011; Betancourt and
Girolami,, 2015; Hoffman et al.,, 2014)-within-Gibbs algorithm, with HMC used
to sample the interaction parameters $\mathbf{\Psi}_{uv}$ for $1\leq u<v\leq
p$, and other parameters updated in Gibbs steps. Although this approach is
highly effective in our experiments, it is important to carefully choose the
step size $e_{0}$ and step length $L_{0}$ in HMC. In practice, we found
$e_{0}\approx 0.01$ and $L_{0}\approx 10$ to work well.

The augmented prior distribution $\pi$ introduced in Section 3.2 creates
difficulty when sampling from the posterior using MCMC. This is because the
normalizing constant

$Z(\mathbf{\Psi},\kappa,\tau_{1},\tau_{2},\nu)=\int\exp\\{-\kappa\widetilde{\mathcal{Q}}(\mathbf{\Psi})\\}\,\pi_{0}(\mathbf{\Psi}\mid\tau_{1},\tau_{2},\nu)\,d\mathbf{\Psi}$

is not easy to evaluate numerically. To bypass this difficulty, we proceed as
in Rao et al., (2016), by augmenting rejected proposals originating from a
rejection sampling mechanism. Introducing the rejected proposals eliminates
the normalizing constant $Z$ and makes sampling from the posterior feasible
using standard MCMC algorithms.

Let $\mathbf{y}=(y_{1},\ldots,y_{n})^{\text{T}}\in\mathbb{R}^{n}$ and
$\mathbf{Z}=[\mathbf{z}_{1}\mid\ldots\mid\mathbf{z}_{n}]^{\text{T}}\in\mathbb{R}^{n\times
q}$. Following Sections 3.2 and 3.4, let
$B_{1},\ldots,B_{p}\in\mathbb{R}^{n\times d}$ and
$S_{1},\ldots,S_{p}\in\mathbb{R}^{n\times m}$ be such that the $i$th row of
$B_{j}$ and $S_{j}$ is given by $\mathbf{b}_{j}^{\text{T}}(x_{ij})$ and
$\mathbf{s}_{j}^{\text{T}}(x_{ij})$, respectively. Let $\mathcal{M}$ be the
block diagonal matrix given by $\mathcal{M}=\mbox{block-
diag}(10^{4},10^{4}\,\mathbb{I}_{q},\Sigma_{M}/\lambda_{1}\ldots,\Sigma_{M}/\lambda_{p})$.
We also let $\mathbf{\widetilde{B}}=[\mathbbm{1}_{n}\mid\mathbf{Z}\mid
B_{1}\mid\ldots\mid B_{p}]\in\mathbb{R}^{n\times(1+q+pd)}$,
$\mathcal{G}=(\alpha,\eta^{\text{T}},\mathbf{\Gamma}_{1}^{\text{T}},\ldots,\mathbf{\Gamma}_{p}^{\text{T}})^{\text{T}}$,
and for $1\leq u<v\leq p$, we let
$\mathbf{h}_{uv}=(S_{u}\theta_{uv,1})^{2}(S_{v}\phi_{uv,1})^{2}-(S_{u}\theta_{uv,2})^{2}(S_{v}\phi_{uv,2})^{2}$
denote the vector of the $uv$th interaction evaluated at the data points, and
$\mathbf{\Theta}_{uv}=(\mathbf{\Psi}_{uv}^{\text{T}},\tau_{uv,1},\tau_{uv,2},\kappa_{uv})^{\text{T}}$
be the set of parameters for the $uv$-th interaction. For a matrix $M_{0}$ and
a finite set $S_{0}$, denote their trace and cardinality by $\mbox{tr}(M_{0})$
and $|S_{0}|$, respectively. Posterior sampling then proceeds as follows.

  1. 1.

Sample $\mathcal{G}\mid-\sim
N\left(\mathcal{A}^{-1}\dfrac{\mathbf{\widetilde{B}}^{\text{T}}\xi}{\sigma^{2}},\mathcal{A}^{-1}\right)$
using Rue, (2001), where

$\mathcal{A}=\dfrac{\mathbf{\widetilde{B}}^{\text{T}}\mathbf{\widetilde{B}}}{\sigma^{2}}\
+\mathcal{M}^{-1},\quad\mathbf{\xi}=\mathbf{y}-\sum_{1\leq u<v\leq
p}\mathbf{h}_{uv}.$

  2. 2.

For each $j=1,\ldots,p$, sample $\lambda_{j}\mid-\sim
G\left(a+\dfrac{d}{2},a+\dfrac{\mathbf{\Gamma}_{j}^{\text{T}}\Sigma_{M}^{-1}\mathbf{\Gamma}_{j}}{2}\right)$.

  3. 3.

For each $u,v$ with $1\leq u<v\leq p$, let
$\mathcal{I}_{uv}=\\{(i_{1},i_{2}):1\leq i_{1}<i_{2}\leq
p,(i_{1},i_{2})\neq(u,v)\\}$.

    1. (a)

Define
$\mathbf{\Delta}_{uv}=\mathbf{y}-\left(\mathbf{\widetilde{B}}\mathcal{G}+\displaystyle\sum_{(u^{\prime},v^{\prime})\in\mathcal{I}_{uv}}\mathbf{h}_{u^{\prime}v^{\prime}}\right)$.

    2. (b)

Given $\kappa_{uv},\tau_{uv,1},\tau_{uv,2},\nu$, repeatedly sample independent
$\mathcal{Y}_{uv,j}\sim\pi_{0}(\cdot\mid\tau_{uv,1},\tau_{uv,2},\nu)$ for
$j=1,2,\ldots$ until a sample $\mathcal{Y}_{uv,j_{0}}$ is accepted; accept
$\mathcal{Y}_{uv,j}$ with probability
$\exp\\{-\kappa_{uv}\widetilde{\mathcal{Q}}(\mathcal{Y}_{uv,j})\\}$ for
$j=1,2,\ldots$. Discard $\mathcal{Y}_{uv,j_{0}}$ and form the set of rejected
proposals
$\mathcal{R}_{uv}=\\{\mathcal{Y}_{uv,1},\ldots,\mathcal{Y}_{uv,|\mathcal{R}_{uv}|}\\}=\\{\mathcal{Y}_{uv,1},\ldots,\mathcal{Y}_{uv,j_{0}-1}\\}$.

    3. (c)

Given $\mathbf{\Delta}_{uv},\mathcal{R}_{uv}$, and $\nu$, use HMC to draw one
sample of $\mathbf{\Theta}_{uv}$, targeting

$\displaystyle\Pi(\mathbf{\Theta}_{uv}\mid\mathbf{\Delta}_{uv},\mathcal{R}_{uv},\nu)\propto$
$\displaystyle\,\pi_{0}(\mathbf{\Psi}_{uv}\mid\tau_{uv,1},\tau_{uv,2},\nu)\,\exp\\{-\kappa_{uv}\widetilde{\mathcal{Q}}(\mathbf{\Psi}_{uv})\\}$
$\displaystyle\prod_{j=1}^{|\mathcal{R}_{uv}|}\pi_{0}(\mathcal{Y}_{uv,j}\mid\tau_{uv,1},\tau_{uv,2},\nu)\left[1-\exp\\{-\kappa_{uv}\widetilde{\mathcal{Q}}(\mathcal{Y}_{uv,j})\\}\right]$
$\displaystyle\pi_{\tau}(\tau_{uv,1})\,\pi_{\tau}(\tau_{uv,2})\,\pi_{\kappa}(\kappa_{uv})\,N(\mathbf{\Delta}_{uv}\mid\mathbf{h}_{uv},\sigma^{2}\mathbb{I}_{n}),$

where $\pi_{\tau}$ is the $C^{+}(0,1)$ density, and $\pi_{\kappa}$ is the
standard lognormal density.

    4. (d)

Update the $uv$-th interaction
$\mathbf{h}_{uv}=(S_{u}\theta_{uv,1})^{2}(S_{v}\phi_{uv,1})^{2}-(S_{u}\theta_{uv,2})^{2}(S_{v}\phi_{uv,2})^{2}$.

  4. 4.

Following Makalic and Schmidt, (2015), we introduce $W$ such that $\nu^{2}\mid
W\sim\mbox{IG}(1/2,1/W)$ and $W\sim\mbox{IG}(1/2,1)$, which leads to $\nu\sim
C^{+}(0,1)$. For each $1\leq u<v\leq p$, we decompose each rejected proposal
$\mathcal{Y}_{uv,j}$ as
$\mathcal{Y}_{uv,j}=(\widetilde{\mathcal{Y}}_{uv,j1}^{\text{T}},\widetilde{\mathcal{Y}}_{uv,j2}^{\text{T}},\widetilde{\mathcal{Y}}_{uv,j3}^{\text{T}},\widetilde{\mathcal{Y}}_{uv,j4}^{\text{T}})^{\text{T}},$
for $j=1,\ldots,|\mathcal{R}_{uv}|$. Form
$\mathbf{R}_{uv,l}\in\mathbb{R}^{|\mathcal{R}_{uv}|\times m}$ with its $j$-th
row given by $\widetilde{\mathcal{Y}}_{uv,jl}^{\text{T}}$, for
$j=1,\ldots,|\mathcal{R}_{uv}|$ and $l=1,2,3,4$. Let $n_{R}=\sum_{1\leq
u<v\leq p}|\mathcal{R}_{uv}|$. The full conditional updates for $\nu^{2}$ and
$W$ are given by:

    1. (a)

$\nu^{2}\mid-\sim\mbox{IG}\left(\dfrac{1}{2}+2m\left\\{\dbinom{p}{2}+n_{R}\right\\},\dfrac{1}{W}+\dfrac{1}{2}\underset{1\leq
u<v\leq p}{\displaystyle\sum}(r_{uv}+t_{uv})\right)$, where

$r_{uv}=\dfrac{\theta_{uv,1}^{\text{T}}\Sigma_{0}^{-1}\theta_{uv,1}+\phi_{uv,1}^{\text{T}}\Sigma_{0}^{-1}\phi_{uv,1}}{\tau_{uv,1}^{2}}+\dfrac{\theta_{uv,2}^{\text{T}}\Sigma_{0}^{-1}\theta_{uv,2}+\phi_{uv,2}^{\text{T}}\Sigma_{0}^{-1}\phi_{uv,2}}{\tau_{uv,2}^{2}},$
$t_{uv}=\dfrac{\mbox{tr}\left\\{\Sigma_{0}^{-1}\left(\mathbf{R}_{uv,1}^{\text{T}}\mathbf{R}_{uv,1}+\mathbf{R}_{uv,2}^{\text{T}}\mathbf{R}_{uv,2}\right)\right\\}}{\tau_{uv,1}^{2}}+\dfrac{\mbox{tr}\left\\{\Sigma_{0}^{-1}\left(\mathbf{R}_{uv,3}^{\text{T}}\mathbf{R}_{uv,3}+\mathbf{R}_{uv,4}^{\text{T}}\mathbf{R}_{uv,4}\right)\right\\}}{\tau_{uv,2}^{2}}.$

    2. (b)

$W\mid-\sim\mbox{IG}(1,1+\nu^{-2}).$

  5. 5.

Sample
$\sigma^{2}\mid-\sim\mbox{IG}\left(\dfrac{n}{2},\dfrac{(\mathbf{y}-\mathbf{\mu})^{\text{T}}(\mathbf{y}-\mathbf{\mu})}{2}\right),$
where $\mu=\mathbf{\widetilde{B}}\mathcal{G}+\displaystyle\sum_{1\leq u<v\leq
p}\mathbf{h}_{uv}.$

We repeat Steps 1-5 for a large number of iterations, and base inference on
the resulting samples after discarding a burn-in. We provide further details
in Section E of the Supplementary Material.

##  4 Simulation Examples

###  4.1 Preliminaries

We carry out simulation studies comparing SAID with the competitors BKMR (Bobb
et al.,, 2015), MixSelect (Ferrari and Dunson,, 2020), HierNet (Bien et al.,,
2013), Family (Haris et al.,, 2016), PIE (Wang et al.,, 2019), and RAMP (Hao
et al.,, 2018). We generate data according to:

$y_{i}=H({\bf x}_{i})+\epsilon_{i},\quad\epsilon_{i}\sim N(0,\sigma_{0}^{2}),$
(14)

where $H$ is the exposure dose response surface decomposed as in (2), and
$\sigma_{0}^{2}$ is the true error variance. As before, we assume that the
exposures satisfy $\mathbf{x}_{i}\in[0,1]^{p}$. We consider two dimensions:
$p=2$ and $p=10$. For $p=2$, we compare the methods in terms of their
estimation performance for varying interaction signal strength and error
variance, across three different types of interactions. For $p=10$, we
consider accuracy in estimation and variable selection for pairwise
interactions.

HierNet, Family, PIE, and RAMP estimate the dose response surface $H$ using
quadratic regression, which also provides estimates of pairwise interaction
functions using bilinear surfaces of the form $f(x,x^{\prime})=\gamma
xx^{\prime}$. MixSelect combines quadratic regression with an additional
nonlinear deviation term in the model, assumed to be orthogonal to the
quadratic regression. This approach provides estimates of pairwise
interactions in the spirit of quadratic regression, while improving
flexibility in capturing $H$. BKMR estimates $H$ using unconstrained Gaussian
processes incorporating variable selection. Thus, we do not consider BKMR as a
competitor when estimating pairwise interactions, as BKMR does not provide
estimates of these interactions.

Given training points $\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n}$, we first
estimate $H$ and the pairwise interactions $h_{uv}$ using the relevant method
and denote the estimates by $\widehat{H}$ and $\widehat{h}_{uv}$,
respectively. For the Bayesian methods, we use the posterior mean as the
estimator. We evaluate the true surface $H$ and true interactions $h_{uv}$ at
test exposure points
$\\{\tilde{\mathbf{x}}_{1},\ldots,\tilde{\mathbf{x}}_{n_{t}}\\}$, denoted by
$\widetilde{\mathbf{H}}=(H(\tilde{\mathbf{x}}_{1}),\ldots,H(\tilde{\mathbf{x}}_{n_{t}}))^{\text{T}}$
and
$\widetilde{\mathbf{h}}_{uv}=(h_{uv}(\tilde{x}_{1u},\tilde{x}_{1v}),\ldots,h_{uv}(\tilde{x}_{n_{t},u},\tilde{x}_{n_{t},v}))^{\text{T}}$,
respectively. Finally, we estimate $\widetilde{\mathbf{H}}$ and
$\widetilde{\mathbf{h}}_{uv}$ for $1\leq u<v\leq p$ by replacing $H$ and
$h_{uv}$ with $\widehat{H}$ and $\widehat{h}_{uv}$, respectively. The
estimates are evaluated using root mean squared error (RMSE).

We also consider the variable selection accuracy of competing methods. For
each method, we estimate the case $1$ and case $2$ errors from classifying an
interaction as synergistic, antagonistic, or null. For any category, the case
$1$ error probability is given by the probability of misclassifying an
interaction as not belonging to that category, when in truth the interaction
is in that category. The case $2$ error probability is given by the
probability of misclassifying an interaction as belonging to that category
when in truth, it does not belong to that category.

###  4.2 Two exposures

We first carry out simulation experiments with $p=2$ exposures. We consider
three different scenarios and evaluate the competitors on test set RMSE in
estimating $H$ and the pairwise interaction $h_{12}$. In all three scenarios,
we train the methods on a sample of size $n=500$ and evaluate on $n_{t}=500$
test points. The pairwise interaction $h_{12}$ is taken to be of the form
$h_{12}(x_{1},x_{2})=\gamma_{0}f(x_{1},x_{2})$, where $\gamma_{0}$ is the
interaction strength and $f$ is a baseline interaction function which we vary
across the three different scenarios. We vary $\gamma_{0}\in\\{1,2\\}$ and
$\sigma_{0}^{2}\in\\{0.1,0.5\\}$, thereby considering low-to-moderate
interaction strength and error variance. For each method in a particular
scenario and a pair $(\gamma_{0},\sigma_{0}^{2})$, we replicate the experiment
$R=10$ times and report the average error across replicates. When fitting
SAID, we assume that the main effects satisfy the origin constraint.

We first consider a scenario where the true interaction is synergistic and
nonlinear in nature and denote this case by SN. To elaborate, the dose
response surface $H(\cdot)$ is given by

$H(x_{1},x_{2})=0.5+x_{1}^{2}+x_{2}^{2}+\gamma_{0}x_{1}^{2}x_{2}^{2}.$

The pairwise interaction $h_{12}(x_{1},x_{2})$ in this case is
$h_{12}(x_{1},x_{2})=\gamma_{0}x_{1}^{2}x_{2}^{2}$. The RMSEs of each method
for varying signal and noise components are provided in Tables 1 and 2.
Although methods such as BKMR and RAMP have similar out-of-sample RMSEs with
respect to SAID, the proposed approach shows superior performance in
estimating the interaction term. In particular, the gains of using SAID are
the most prominent when the signal strength of the interaction is the highest,
indicating lack of flexibility of the other methods.

Signal and Noise | SAID | BKMR | MixSelect | HierNet | Family | PIE | RAMP
---|---|---|---|---|---|---|---
$\gamma_{0}=1,\sigma_{0}^{2}=0.1$ | 0.05 | 0.04 | 0.15 | 0.22 | 0.40 | 1.41 | 0.05
$\gamma_{0}=1,\sigma_{0}^{2}=0.5$ | 0.10 | 0.10 | 0.16 | 0.22 | 0.41 | 1.29 | 0.09
$\gamma_{0}=2,\sigma_{0}^{2}=0.1$ | 0.05 | 0.05 | 0.13 | 0.24 | 0.51 | 2.18 | 0.07
$\gamma_{0}=2,\sigma_{0}^{2}=0.5$ | 0.10 | 0.10 | 0.19 | 0.24 | 0.51 | 2.16 | 0.10
Table 1: RMSE of competitors when estimating H in case SN. Here $n=500$ and $p=2$. Signal and Noise | SAID | MixSelect | HierNet | Family | PIE | RAMP
---|---|---|---|---|---|---
$\gamma_{0}=1,\sigma_{0}^{2}=0.1$ | 0.06 | 0.16 | 0.19 | 0.18 | 0.18 | 0.18
$\gamma_{0}=1,\sigma_{0}^{2}=0.5$ | 0.11 | 0.19 | 0.16 | 0.18 | 0.24 | 0.21
$\gamma_{0}=2,\sigma_{0}^{2}=0.1$ | 0.08 | 0.24 | 0.17 | 0.37 | 0.34 | 0.34
$\gamma_{0}=2,\sigma_{0}^{2}=0.5$ | 0.13 | 0.28 | 0.18 | 0.37 | 0.37 | 0.36
Table 2: RMSE of competitors when estimating $h_{12}$ in case SN. Here $n=500$
and $p=2$.

Secondly, we assume the data generating process has linear main effects and a
synergistic linear interaction, namely a quadratic regression (QR) setup. We
let $H$ be

$H(x_{1},x_{2})=0.5+x_{1}+x_{2}+\gamma_{0}x_{1}x_{2}.$

The interaction is given by $h_{12}(x_{1},x_{2})=\gamma_{0}x_{1}x_{2}$.
Lastly, we consider a scenario where the interaction is nonlinear and is
neither synergistic, antagonistic, or null, while keeping the main effects as
before. We call this the mis-specified interaction (MIS) case and let $H$ be

$H(x_{1},x_{2})=0.5+x_{1}+x_{2}+\gamma_{0}(x_{1}x_{2}-2x_{1}^{2}x_{2}^{2}).$

In this case, the interaction is given by
$h_{12}(x_{1},x_{2})=\gamma_{0}(x_{1}x_{2}-2x_{1}^{2}x_{2}^{2}).$ The
interaction $h_{12}$ is non-negative for $x_{1}x_{2}\leq 1/2$ and non-positive
for $x_{1}x_{2}\geq 1/2$, and thus is neither synergistic, antagonistic, or
null. To ease exposition, we defer the results obtained from scenarios QR and
MIS to Section G of the Supplementary Material. In the scenario QR, quadratic
regression approaches such as RAMP and MixSelect perform better as the data is
also generated from a quadratic regression setup. However, the performance of
SAID is similar to its competitors for estimating both the dose response
surface and the interaction surface. In the scenario MIS, BKMR and SAID
perform the best in terms of estimating $H$. When estimating the interaction,
SAID outperforms its competitors due to the lack of flexibility for other
pairwise interaction approaches.

###  4.3 More than two exposures

In this subsection, we consider a moderate dimensional example with $p=10$ and
thus $\binom{10}{2}=45$ pairwise interactions. We assume that the true data
generating model has $5$ synergistic, $5$ antagonistic, and $35$ null pairwise
interactions. We let $\sigma_{0}^{2}=0.2$ with the number of training and test
points taken to be $n=1000$ and $n_{t}=500$, respectively. The experiment is
replicated $R=20$ times. We fit all the quadratic regression methods assuming
weak heredity of the pairwise interactions. Family is not considered due to
unstable estimates; furthermore, except when estimating $H$, we do not
consider BKMR.

We assume (14) and decompose $H({\bf x})$ in the following way:

$H({\bf x})=\alpha_{0}+M_{0}({\bf x})+S_{0}({\bf x})+A_{0}({\bf x}),$

where

$\displaystyle\alpha_{0}$ $\displaystyle=-5/6$ $\displaystyle M_{0}({\bf x})$
$\displaystyle=\left(x_{1}+x_{1}^{2}\right)+\dfrac{x_{2}}{2}+x_{7}^{3},$
$\displaystyle S_{0}({\bf x})$
$\displaystyle=4(x_{1}-x_{1}^{2})x_{2}+x_{1}x_{9}+x_{2}^{2}x_{3}^{2}+x_{3}x_{8}+\dfrac{(e^{x_{5}}-1)x_{10}}{e-1},$
$\displaystyle A_{0}({\bf x})$
$\displaystyle=-\left[x_{1}x_{3}+x_{2}^{2}x_{5}+\dfrac{27}{4}x_{4}^{2}(1-x_{4})x_{9}+x_{7}x_{10}+x_{8}x_{9}^{2}\right].$

Here, $M_{0}$, $S_{0}$, $A_{0}$ denote the true main effects, synergistic
interaction effects, and the antagonistic interaction effects, respectively.
The form of the interactions include pairwise linear, nonlinear polynomial,
and nonlinear interactions which are not of polynomial form. Each pairwise
interaction has absolute maximum value equal to $1$. As before, we assume the
main effects start from the origin when fitting SAID. The results comparing
the methods are given in Table 4. For comparison, the RMSE when estimating $H$
using BKMR is 0.17.

In terms of estimation accuracy, SAID performs uniformly better than its
competitors in estimating both the overall surface $H$ and the interaction
surface $I$, given by $I({\bf x})=S_{0}({\bf x})+A_{0}({\bf x}).$ We believe
that this is due to the SAID prior on the pairwise interactions being flexible
enough to capture nonlinear interactions while also efficiently extracting
interaction signal in the presence of noise. Furthermore, the SAID framework
does not require any heredity assumptions, helping detection of pairwise
interactions even in the absence of main effects of one or both of the
corresponding exposures. In terms of variable selection, the case 1 error
probabilities for synergistic and antagonistic detections are similar across
the methods. However, SAID has considerably less case $1$ error when detecting
null interactions, indicating that SAID more accurately classifies null
interactions to be null compared with its competitors. SAID also shows
superior performance in terms of case $2$ error when classifying both
synergistic and antagonistic interactions. This indicates that SAID has a
lower probability of misclassifying an interaction as synergistic or
antagonistic when it is not so, compared with its competitors. Out of the
methods considered, SAID is the only one with case 1 and case 2 classification
errors less than $0.05$ for each case. In terms of uncertainty quantification
of interactions evaluated at the test points, the $95\%$ posterior credible
intervals obtained from SAID have $94.6\%$ coverage, averaged over $R=20$
replicates and all $45$ interactions.

| SAID | MixSelect | HierNet | PIE | RAMP
---|---|---|---|---|---
Overall Surface RMSE | 0.13 | 0.29 | 0.30 | 0.75 | 0.18
Interaction Surface RMSE | 0.21 | 0.52 | 0.48 | 0.42 | 0.43
Synergistic Case 1 Probability | 0.008 | 0.040 | 0.001 | 0.008 | 0.015
Synergistic Case2 Probability | 0.001 | 0.220 | 0.450 | 0.210 | 0.230
Antagonistic Case1 Probability | 0.005 | 0.013 | 0.002 | 0.006 | 0.008
Antagonistic Case2 Probability | 0.001 | 0.050 | 0.240 | 0.040 | 0.050
Null Case1 Probability | 0.001 | 0.135 | 0.345 | 0.125 | 0.140
Null Case2 Probability | 0.014 | 0.056 | 0.002 | 0.016 | 0.026
Table 3: Comparison of the methods in terms of RMSE and variable selection
accuracy for $p=10$.

##  5 Analysis of Kidney Function Data

###  5.1 Preliminaries

In this Section, we apply the proposed Synergistic Antagonistic Interaction
Detection (SAID) approach to the NHANES 2015-16 data. We are interested in
detecting synergistic and antagonistic interactions between heavy metals
affecting kidney function. Following Section 2 and the results in Section 4,
it is evident that the current methods either do not provide inferences on
interactions or lack flexibility in characterizing interactions. We now
illustrate how SAID detects synergistic, antagonistic, and null interactions.

As discussed in Section 2.1, we assess kidney function of individuals through
their urine creatinine (uCr) levels, measured in mg/dL. Heavy metal
concentrations are also measured in urine. We consider $13$ heavy metals,
namely Antimony (Sb), Barium (Ba), Cadmium (Cd), Cesium (Cs), Cobalt (Co),
Lead (Pb), Manganese (Mn), Molybdenum (Mo), Strontium (Sr), Thallium (Tl), Tin
(Sn), Tungsten (W), and Uranium (U), all measured in $\mu$g/mL. Following
Section 2.2, we remove individuals with albuminuria and missing entries from
the original data set. The sample size after removal of such entries is
$n=1979$. As described in Section 2.3, we multiply both the uCr and the heavy
metal concentration levels by the individual-specific urine flow rate to
obtain their dilution-adjusted versions. Furthermore, we also consider age (in
years), sex (0 for males and 1 for females), ethnicity (Non-Hispanic White,
Non-Hispanic Black, Mexican American, Other Hispanic, and Other), and body
mass index (BMI) of the subject as covariates.

Before beginning our analysis, we first (natural) log-transformed the
dilution-adjusted urine creatinine levels. We use a marginal cumulative
distribution function (CDF) transformation for standardizing exposure
variables. Suppose the dilution-adjusted exposure variables are denoted by
$E_{1},\ldots,E_{p}$, with the marginal CDF of variable $j$ denoted by $F_{j}$
for $j=1,\ldots,p$. To estimate $F_{j}$, we use kernel density estimation to
first estimate the marginal density of $E_{j}$ and then estimate the induced
CDF $\hat{F}_{j}$ from the density estimate. The transformed exposures are
defined as $x_{j}=\hat{F}_{j}(E_{j})\in[0,1]$. Standardizing exposures in this
manner facilitates statistical inferences by reducing the tendency for the
exposure data to be unevenly distributed, with very sparse observations for
certain ranges of exposure. Transforming the exposures does not complicate
interpretation, as we can convert dose response surfaces back to the original
units. In addition, the transformed exposures are directly interpretable as
quantiles of the exposure distribution in the sample; for example, $x_{j}=0.5$
reflects a median value for the $j$th exposure.

For the $i$-th individual in the data set, let $y_{i}\in\mathbb{R}$ be the log
of the dilution-adjusted urinary creatinine level,
$\mathbf{x}_{i}\in[0,1]^{13}$ be the $13$ dilution-adjusted urinary metal
concentrations after CDF transformations, and
$\mathbf{z}_{i}\in\mathbb{R}^{7}$ be the covariate vector, with dummy
variables for the different ethnic groups. We use “Non-Hispanic White” as the
baseline category in defining indicators. The dimensions of exposures and
covariates are $p=13$ and $q=7$, respectively. We also standardized the
covariates age and BMI to have variance $1$ before fitting the model.

Following Section 3.5, we employed a Hamiltonian Monte Carlo (HMC)-within-
Gibbs sampler to obtain posterior samples of the model parameters. We ran the
sampler for a total of $15000$ iterations and discarded the first $5000$
iterations as burn-in. We standardized the $y$ to have variance $1$ before
fitting the model and then used an integral cutoff of $0.01$ as described in
Section 3.3 in order to compute PIPs of pairwise interactions. After computing
the PIPs, we present our inferences in the original scale, that is, where $y$
is not standardized. This is achieved simply by multiplying the estimates of
the intercept, main effects, interaction effects, covariate effects, and
measurement error standard deviation obtained from the fitted model by the
standard deviation of $y$.

###  5.2 Model Diagnostics

We first assess MCMC convergence. To improve mixing, we perturb the step-size
$e_{0}$ by a small factor every $500$ iterations. As a measure of mixing of
the chain, we look at the MCMC samples of the error variance $\sigma^{2}$. The
effective sample size of $\sigma^{2}$ is $33.7\%$ of the total number of MCMC
samples, after discarding the burn-in. We found this proportion to remain
fairly robust across longer chains. On a MacBook Pro with M1 CPU and 16 GB of
RAM, it took $\sim 40$ minutes to fit SAID on this dataset; for comparison,
MixSelect took $\sim 6$ hours with the same number of MCMC iterates.
Computation time can be improved by taking shorter chains, as $15000$ was
conservative based on our assessments.

Next, we assess goodness-of-fit by inspecting standardized residuals and
carrying out posterior predictive checks (Gelman et al.,, 1996). We generate
posterior predictive samples corresponding to each training point. To validate
the assumption of normality of the measurement error, we look at a Q-Q plot of
standardized residuals, defined to be the difference of the observed and
predicted responses and then standardized. We also compare the marginal
density of the observed response variable with the marginal densities of $3$
randomly chosen MCMC samples of predicted responses. We provide both the Q-Q
plot of the obtained standardized residuals and the comparison of marginal
density plots in Figure 3. Figure 3(a) indicates that the normality assumption
is justified; furthermore, the densities of the randomly chosen predicted
response samples in Figure 3(b) very closely resemble the marginal density of
the observed response. Lastly, the coverage of $95\%$ posterior predictive
intervals of the responses, averaged over all observed responses, is $96\%$.

|
---|---
(a) Q-Q Plot of standardized residuals. | (b) Marginal density plot of generated responses.
Figure 3: Figure showing Q-Q plot of standardized residuals and marginal
density plots of posterior predictive samples of urine creatinine levels,
obtained from the NHANES 2015-16 data. Deviation from the reference blue line
in Q-Q plot is deviation from normality.

###  5.3 Results

In this subsection, we discuss the main results of our analyses of the data
described in Section 5.1. For each exposure, we constrain the main effects to
start from $0$ at the minimum observed value of that exposure. The $95\%$
posterior credible interval of the error variance $\sigma^{2}$ is
$[0.080,0.091]$ with a posterior mean of $0.086$. The exposures and covariates
together explain $89\%$ of the variation in the response.

We observed nonlinear main effects for Antimony, Cadmium, Cobalt, Cesium,
Molybdenum, Tin, and Uranium. The general trend for the main effects of these
heavy metals is an increase in log dilution-adjusted urine creatinine as
exposure to metal concentrations increases. Since urine concentrations of
creatinine are directly correlated with serum concentrations of creatinine,
this might indicate higher kidney stress at higher levels of metal exposure.
The main effects have either a monotonic or hill-shaped pattern, which is as
expected in studying health effects of potentially toxic exposures. We found
exposure to high doses of Cesium to increase log dilution-adjusted urine
creatinine the most, followed by Cadmium. For illustration, we plot the main
effects of Cadmium, Cesium, Molybdenum, and Uranium in Figure 4, with both the
response and the exposures shown in log scale. We provide further plots for
the main effects of the other exposures in Section H of the Supplementary
Material. When the exposures are CDF transformed, we can evaluate the main
effect of an exposure at a desired quantile. As an illustration, the metals
Antimony, Cadmium, Cobalt, Cesium, Molybdenum, Tin, and Uranium at their
median exposure levels have a main effect of $0.23,0.69,0.52,0.88,0.39,0.22$,
and $0.22$, respectively, on log dilution-adjusted Creatinine. Similar main
effects of heavy metal exposures on kidney function have been detected in
earlier literature. For example, Ferraro et al., (2010) found Cadmium exposure
to be associated with increased risk of chronic kidney disease, based on an
analysis of NHANES data from 1999-2006. In a recent study, Rahman et al.,
(2022) found Cesium, Cadmium, and Antimony to be associated with kidney
damage.

Figure 4: Plots showing the main effects of dilution-adjusted Cadmium, Cesium,
Molybdenum, and Uranium on log dilution-adjusted Creatinine. Exposure levels
are in log scale. Black line denotes posterior mean and shaded regions denote
pointwise $95\%$ posterior credible intervals.

The proposed approach detects multiple synergistic and antagonistic
interactions between exposures. Excluding the interactions with
$\mbox{PIP}<0.5$, the detected interactions were Cadmium and Tin (PIP $>$
0.99), Molybdenum and Tin (PIP $>$ 0.99), Cadmium and Manganese (PIP = 0.99),
and Cobalt and Manganese (PIP = 0.98). For these interactions, we also compute
the posterior synergistic probability (PSP) and posterior antagonistic
probability (PAP). The interactions Cadmium and Tin (PSP $>$ 0.99) and Cobalt
and Manganese (PSP = 0.98) have high posterior probability of being
synergistic. On the other hand, the interactions Molybdenum and Tin (PAP $>$
0.99) and Cadmium and Manganese (PAP = 0.99) have high posterior probabilities
of being antagonistic. The interactions typically demonstrate flat behavior
for most of the exposure domain and synergy/antagonism for the rest of the
domain. Nonlinear surfaces of this kind cannot be captured by quadratic
regression approaches. We believe this to be a possible reason behind
MixSelect being unable to detect these interactions. In Figure 5, we plot the
interaction surfaces of Cadmium and Tin, Cobalt and Manganese, and Cadmium and
Manganese. The plots for the interaction between Molybdenum and Tin can be
found in Section H of the Supplementary Material. Each row in Figure 5, from
left to right, shows the pointwise $2.5\%$ quantile, mean, and pointwise
$97.5\%$ quantile of the posterior samples of the pairwise interaction,
evaluated on a $30\times 30$ regular grid of points across the exposure
values.

Regarding the covariate effects, we detected a negative association with age
with an estimated coefficient $-0.10$ and $95\%$ CI given by $[-0.12,-0.08]$.
As expected, females had significantly lower urine creatinine levels than
males, with the estimated coefficient of sex being $-0.25$ with $95\%$ CI
given by $[-0.28,-0.22]$. In addition, average creatinine levels increased
with BMI, with the coefficient on BMI estimated as $0.11$ with $95\%$ CI of
$[0.10,0.13]$. The ethnicity category “Non-Hispanic Black” was found to have
higher log dilution-adjusted urine creatinine concentrations than that of the
baseline category “Non-Hispanic White”, with an estimated coefficient of
$0.17$ and a $95\%$ CI $[0.13,0.21]$. The ethnicity categories “Mexican
American”, “Other Hispanic”, and “Other” were found to have lower log
dilution-adjusted urine creatinine concentrations than that of “Non-Hispanic
White”; their estimated coefficients are $-0.12$ with $95\%$ CI given by
$[-0.16,-0.08]$, $-0.05$ with $95\%$ CI given by $[-0.10,-0.01]$, and $-0.11$
with a $95\%$ CI given by $[-0.16,-0.07]$, respectively. Similar observations
have been made previously in literature; refer to James et al., (1988).

Figure 5: Plots showing the interaction effects of dilution-adjusted Cadmium
and Tin, Cobalt and Manganese, and Cadmium and Manganese on dilution-adjusted
log Creatinine. Exposure levels are in log scale. Each plot shows the
pointwise $2.5\%$ posterior credible surface, the posterior mean, and the
$97.5\%$ posterior credible surface from left to right.

##  6 Discussion

In this paper, we proposed a framework for inference on main effects and
pairwise synergistic, antagonistic, or null interaction effects, motivated by
the mixtures problem in environmental epidemiology. Instead of sharply
imposing constraints, we use a shrinkage prior that penalizes deviations from
synergistic, antagonistic, or null interactions. When the interaction is
synergistic or antagonistic, the proposed approach identifies the correct
direction of interaction, while greatly relaxing the penalty when the
interaction is neither synergistic nor antagonistic. We proposed a variable
selection approach to identify non-null interactions and additionally provide
posterior probabilities of an interaction being synergistic and antagonistic.
Unlike most quadratic regression approaches, we do not rely on heredity
assumptions to detect interactions. The proposed approach improves upon
existing methods by providing greater flexibility over quadratic regression,
while improving statistical power and interpretability in detecting
interactions over unrestricted nonparametric models. We demonstrated these
improvements over existing state-of-the-art methods both from the perspective
of estimation and variable selection.

In our NHANES analysis, we found a variety of significant interactions among
metal exposures impacting kidney function. To our knowledge, these
interactions are currently unknown to the epidemiology and public health
communities. It will be very interesting to validate these results in other
cohorts, to study the mechanisms by which synergistic and antagonistic
interactions occur in these particular metals, and also assess the regulatory
implications. Our approach for assessing interactions in observational
epidemiology data can be used to identifying specific pairs of chemicals to
study in more detail in in vivo and in vitro assays.

Although we assumed the response to be continuous, it is straightforward to
incorporate other types of responses in our framework, such as binary or count
variables. In particular, we can model the response variables as falling in an
exponential family with the linear predictor having exactly the form described
earlier in this article. For binary outcomes and assuming a logistic link
function, our proposed computational algorithm can be easily modified by
relying on Polya-Gamma data augmentation (Polson et al.,, 2013). A related
approach can be used for counts, or our Gibbs steps can be replaced with
gradient-based updates to avoid data augmentation. It is additionally
straightforward to include further applied complications, such as spatially or
temporally dependent data, by incorporating appropriate terms in the linear
predictor.

The exposures in the NHANES 2015-16 data considered in this paper are mild-to-
moderately correlated. In scenarios where the exposures are more heavily
correlated, individually interpreting main effects and interaction effects is
often infeasible. In this scenario, one could think of dividing the exposure
variables into groups based on chemical class, clustering algorithms, or
latent factor models. It is then of interest to investigate presence of
interactions between chemical groups.

##  7 Acknowledgements

R package SAID was used in the application and numerical experiments. This
research was supported by grants R01-ES028804 and R01-ES027498 from the
National Institute of Environmental Health Sciences of the National Institutes
of Health. The authors would like to thank Federico Ferrari for helpful
comments.

## References

  * Barr et al., (2005) Barr, D. B., Wilder, L. C., Caudill, S. P., Gonzalez, A. J., Needham, L. L., and Pirkle, J. L. (2005).  Urinary creatinine concentrations in the us population: implications for urinary biologic monitoring measurements.  Environmental Health Perspectives, 113(2):192–200.
  * Baxmann et al., (2008) Baxmann, A. C., Ahmed, M. S., Marques, N. C., Menon, V. B., Pereira, A. B., Kirsztajn, G. M., and Heilberg, I. P. (2008).  Influence of muscle mass and physical activity on serum and urinary creatinine and serum cystatin c.  Clinical Journal of the American Society of Nephrology, 3(2):348–354.
  * Betancourt and Girolami, (2015) Betancourt, M. and Girolami, M. (2015).  Hamiltonian monte carlo for hierarchical models.  Current Trends in Bayesian Methodology with Applications, 79(30):2–4.
  * Bien et al., (2013) Bien, J., Taylor, J., and Tibshirani, R. (2013).  A lasso for hierarchical interactions.  Annals of Statistics, 41(3):1111.
  * Bobb et al., (2015) Bobb, J. F., Valeri, L., Claus Henn, B., Christiani, D. C., Wright, R. O., Mazumdar, M., Godleski, J. J., and Coull, B. A. (2015).  Bayesian kernel machine regression for estimating the health effects of multi-pollutant mixtures.  Biostatistics, 16(3):493–508.
  * Brezger and Lang, (2006) Brezger, A. and Lang, S. (2006).  Generalized structured additive regression based on bayesian p-splines.  Computational Statistics & Data Analysis, 50(4):967–991.
  * Carvalho et al., (2009) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2009).  Handling sparsity via the horseshoe.  In Artificial Intelligence and Statistics, pages 73–80. PMLR.
  * Chipman et al., (2010) Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).  Bart: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1):266–298.
  * Czarnota et al., (2015) Czarnota, J., Gennings, C., and Wheeler, D. C. (2015).  Assessment of weighted quantile sum regression for modeling chemical mixtures and cancer risk.  Cancer Informatics, 14:CIN–S17295.
  * De Boor, (1978) De Boor, C. (1978).  A Practical Guide to Splines, volume 27.  Springer-Verlag New York.
  * Ferrari and Dunson, (2020) Ferrari, F. and Dunson, D. B. (2020).  Identifying main effects and interactions among exposures using gaussian processes.  The Annals of Applied Statistics, 14(4):1743.
  * Ferrari and Dunson, (2021) Ferrari, F. and Dunson, D. B. (2021).  Bayesian factor analysis for inference on interactions.  Journal of the American Statistical Association, 116(535):1521–1532.
  * Ferraro et al., (2010) Ferraro, P. M., Costanzi, S., Naticchia, A., Sturniolo, A., and Gambaro, G. (2010).  Low level exposure to cadmium increases the risk of chronic kidney disease: analysis of the nhanes 1999-2006.  BMC Public Health, 10(1):1–8.
  * Forbes and Bruining, (1976) Forbes, G. and Bruining, G. J. (1976).  Urinary creatinine excretion and lean body mass.  The American Journal of Clinical Nutrition, 29(12):1359–1366.
  * Gelman et al., (1996) Gelman, A., Meng, X.-L., and Stern, H. (1996).  Posterior predictive assessment of model fitness via realized discrepancies.  Statistica Sinica, pages 733–760.
  * Hao et al., (2018) Hao, N., Feng, Y., and Zhang, H. H. (2018).  Model selection for high-dimensional quadratic regression via regularization.  Journal of the American Statistical Association, 113(522):615–625.
  * Haris et al., (2016) Haris, A., Witten, D., and Simon, N. (2016).  Convex modeling of interactions with strong heredity.  Journal of Computational and Graphical Statistics, 25(4):981–1004.
  * Hays et al., (2015) Hays, S. M., Aylward, L. L., and Blount, B. C. (2015).  Variation in urinary flow rates according to demographic characteristics and body mass index in nhanes: potential confounding of associations between health outcomes and urinary biomarker concentrations.  Environmental Health Perspectives, 123(4):293–300.
  * Hoffman et al., (2014) Hoffman, M. D., Gelman, A., et al. (2014).  The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo.  J. Mach. Learn. Res., 15(1):1593–1623.
  * Jain, (2016) Jain, R. (2016).  Associated complex of urine creatinine, serum creatinine, and chronic kidney disease.  Epidemiology (Sunnyvale), 6(234):2161–1165.
  * James et al., (1988) James, G. D., Sealey, J. E., Alderman, M., Ljungman, S., Mueller, F. B., Pecker, M. S., and Laragh, J. H. (1988).  A longitudinal study of urinary creatinine and creatinine clearance in normal subjects: race, sex, and age differences.  American Journal of Hypertension, 1(2):124–131.
  * Jeng et al., (2021) Jeng, P.-H., Huang, T.-R., Wang, C.-C., and Chen, W.-L. (2021).  Clinical relevance of urine flow rate and exposure to polycyclic aromatic hydrocarbons.  International Journal of Environmental Research and Public Health, 18(10):5372.
  * Joubert et al., (2022) Joubert, B. R., Kioumourtzoglou, M.-A., Chamberlain, T., Chen, H. Y., Gennings, C., Turyk, M. E., Miranda, M. L., Webster, T. F., Ensor, K. B., Dunson, D. B., et al. (2022).  Powering research through innovative methods for mixtures in epidemiology (prime) program: Novel and expanded statistical methods.  International Journal of Environmental Research and Public Health, 19(3):1378.
  * Kashani et al., (2020) Kashani, K., Rosner, M. H., and Ostermann, M. (2020).  Creatinine: from physiology to clinical application.  European Journal of Internal Medicine, 72:9–14.
  * Kass and Raftery, (1995) Kass, R. E. and Raftery, A. E. (1995).  Bayes factors.  Journal of the American Statistical Association, 90(430):773–795.
  * Kim et al., (2015) Kim, N. H., Hyun, Y. Y., Lee, K.-B., Chang, Y., Rhu, S., Oh, K.-H., and Ahn, C. (2015).  Environmental heavy metal exposure and chronic kidney disease in the general population.  Journal of Korean Medical Science, 30(3):272–277.
  * Lang and Brezger, (2004) Lang, S. and Brezger, A. (2004).  Bayesian p-splines.  Journal of Computational and Graphical Statistics, 13(1):183–212.
  * Luo and Hendryx, (2020) Luo, J. and Hendryx, M. (2020).  Metal mixtures and kidney function: An application of machine learning to nhanes data.  Environmental Research, 191:110126.
  * Makalic and Schmidt, (2015) Makalic, E. and Schmidt, D. F. (2015).  A simple sampler for the horseshoe estimator.  IEEE Signal Processing Letters, 23(1):179–182.
  * Middleton et al., (2016) Middleton, D. R., Watts, M. J., Lark, R. M., Milne, C. J., and Polya, D. A. (2016).  Assessing urinary flow rate, creatinine, osmolality and other hydration adjustment methods for urinary biomonitoring using nhanes arsenic, iodine, lead and cadmium data.  Environmental Health, 15(1):1–13.
  * Molitor et al., (2010) Molitor, J., Papathomas, M., Jerrett, M., and Richardson, S. (2010).  Bayesian profile regression with an application to the national survey of children’s health.  Biostatistics, 11(3):484–498.
  * Neal et al., (2011) Neal, R. M. et al. (2011).  Mcmc using hamiltonian dynamics.  Handbook of Markov Chain Monte Carlo, 2(11):2.
  * Pollack et al., (2015) Pollack, A. Z., Mumford, S. L., Mendola, P., Perkins, N. J., Rotman, Y., Wactawski-Wende, J., and Schisterman, E. F. (2015).  Kidney biomarkers associated with blood lead, mercury, and cadmium in premenopausal women: a prospective cohort study.  Journal of Toxicology and Environmental Health, Part A, 78(2):119–131.
  * Polson et al., (2013) Polson, N. G., Scott, J. G., and Windle, J. (2013).  Bayesian inference for logistic models using pólya–gamma latent variables.  Journal of the American Statistical Association, 108(504):1339–1349.
  * Rahman et al., (2022) Rahman, H. H., Niemann, D., and Munson-McGee, S. H. (2022).  Association of albumin to creatinine ratio with urinary arsenic and metal exposure: evidence from nhanes 2015–2016.  International Urology and Nephrology, 54(6):1343–1353.
  * Ramsay, (1988) Ramsay, J. O. (1988).  Monotone regression splines in action.  Statistical Science, 3(4):425 – 441.
  * Rao et al., (2016) Rao, V., Lin, L., and Dunson, D. B. (2016).  Data augmentation for models based on rejection sampling.  Biometrika, 103(2):319–335.
  * Rue, (2001) Rue, H. (2001).  Fast sampling of gaussian markov random fields.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(2):325–338.
  * Shively et al., (2009) Shively, T. S., Sager, T. W., and Walker, S. G. (2009).  A bayesian approach to non-parametric monotone function estimation.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(1):159–175.
  * Wang et al., (2019) Wang, C., Jiang, B., and Zhu, L. (2019).  Penalized interaction estimation for ultrahigh dimensional quadratic regression.  arXiv preprint arXiv:1901.07147.
  * Wei et al., (2020) Wei, R., Reich, B. J., Hoppin, J. A., and Ghosal, S. (2020).  Sparse bayesian additive nonparametric regression with application to health effects of pesticides mixtures.  Statistica Sinica, 30(1):55–79.
  * Williams and Rasmussen, (2006) Williams, C. K. and Rasmussen, C. E. (2006).  Gaussian Processes for Machine Learning, volume 2.  MIT press Cambridge, MA.

## Supplementary Material

##  Appendix A B-spline Functions

For an extensive overview of B-spline basis functions, we refer the reader to
De Boor, (1978). As a default, we choose cubic splines for modeling both the
main and interaction effects throughout the paper; however, one could vary the
degree of the B-splines as required. We now describe the choice of B-splines
for modeling the interaction effects and main effects in Sections 3.2 and 3.4,
respectively.

For the interaction function $h_{uv}$, we let $h_{uv}$ satisfy
$h_{uv}(x_{u},0)=0$ for all $x_{u}\in[0,1]$ and $h_{uv}(0,x_{v})=0$ for all
$x_{v}\in[0,1]$. To achieve this, we first construct the set of B-spline
functions $s_{u0}(x_{u}),\ldots,s_{um}(x_{u})$ along the $u$-th dimension for
$u=1,\ldots,p$, where $s_{u0}$ is the intercept spline. The intercept spline
$s_{u0}$ the only B-spline basis function satisfying $s_{u0}(0)\neq 0$. Thus,
to model the interaction effects, we simply ignore the intercept spline
$s_{u0}$ and only choose the B-spline functions $s_{uj}(\cdot)$ that satisfy
$s_{uj}(0)=0$.

In order to model the main effects, we proceed as follows. Suppose we have
B-spline basis functions $t_{u0},\ldots,t_{ud}$, with $t_{u0}$ being the
intercept spline. If the main effects satisfy the origin constraint, we simply
ignore $t_{u0}$ and let the basis functions be $b_{u,j}(x_{u})=t_{uj}(x_{u})$
for $u=1,\ldots,p$ and $j=1,\ldots,d$. If the main effects satisfy the
integral constraint, we let
$b_{u,j}(x_{u})=t_{uj}(x_{u})-\int_{0}^{1}t_{uj}(x)\,dx$ for $u=1,\ldots,p$
and $j=0,\ldots,d$.

##  Appendix B P-spline Priors

In order to ensure a function $g(x)=\sum_{j=1}^{m}s_{j}(x)\beta_{j}$ expressed
as a linear combination of the B-splines $s_{1},\ldots,s_{m}$ is sufficiently
smooth, we follow the P-spline approach described in Lang and Brezger, (2004).
P-spline priors aim to promote smoothness by shrinking coefficients of
adjacent splines towards each other.

We use a P-spline prior of order $1$ along with a slight modification to make
the it a proper prior distribution. The order $1$ prior described in Lang and
Brezger, (2004) is obtained by considering the hierarchical model
$\beta_{j}\mid\beta_{j-1}\sim N(\beta_{j-1},\tau^{2})$ for $j=2,\ldots,m$, and
then letting $\beta_{1}$ have an improper prior, given by
$\pi(\beta_{1})\propto 1$. However, we avoid the use of improper priors, and
let $\beta_{1}\sim N(0,\tau^{2})$ as well. Under this modified prior
specification, the joint distribution of
$\mathbf{\beta}=(\beta_{1},\ldots,\beta_{m})^{\text{T}}$ is given by

$\pi(\mathbf{\beta})\propto\exp\left[-\dfrac{1}{2\tau^{2}}\left\\{\beta_{1}^{2}+\sum_{j=2}^{m}(\beta_{j}-\beta_{j-1})^{2}\right\\}\right].$

This can be rewritten as $\mathbf{\beta}\sim N(0,\tau^{2}\Sigma_{0})$ for a
positive-definite matrix $\Sigma_{0}$ whose entries are known.

##  Appendix C Model Identifiability

We first state a result which ensures identifiability of the main and
interaction effects in the assumed model.

######  Proposition C.1.

Suppose the dose response function is given by

$H(\mathbf{x})=\alpha+\sum_{j=1}^{p}f_{j}(x_{j})+\sum_{1\leq u<v\leq
p}h_{uv}(x_{u},x_{v}),$

for $\mathbf{x}=(x_{1},\ldots,x_{p})^{\text{T}}\in[0,1]^{p}$, $p\geq 2$.
Assume the following conditions hold:

  1. 1.

$f_{1},\ldots,f_{p}$ together satisfy either (a) or (b), where

    1. (a)

Origin Constraint: $f_{j}(0)=0$ for $j=1,\ldots,p$, or

    2. (b)

Integral Constraint: $\displaystyle\int_{0}^{1}f_{j}(x_{j})\,dx_{j}=0$ for
$j=1,\ldots,p$.

  2. 2.

For all $1\leq u<v\leq p$, $h_{uv}(x_{u},0)=0$ and $h_{uv}(0,x_{v})=0$ for all
$x_{u},x_{v}\in[0,1]$.

If two tuples $\left(\alpha,\\{f_{j}\\}_{j=1}^{p},\\{h_{uv}\\}_{1\leq u<v\leq
p}\right)$ and
$\left(\tilde{\alpha},\\{\tilde{f}_{j}\\}_{j=1}^{p},\\{\tilde{h}_{uv}\\}_{1\leq
u<v\leq p}\right)$ satify

$\alpha+\sum_{j=1}^{p}f_{j}(x_{j})+\sum_{1\leq u<v\leq
p}h_{uv}(x_{u},x_{v})=\tilde{\alpha}+\sum_{j=1}^{p}\tilde{f}_{j}(x_{j})+\sum_{1\leq
u<v\leq p}\tilde{h}_{uv}(x_{u},x_{v}),$

for all $\mathbf{x}\in[0,1]^{p}$, then
$\left(\alpha,\\{f_{j}\\}_{j=1}^{p},\\{h_{uv}\\}_{1\leq u<v\leq
p}\right)=\left(\tilde{\alpha},\\{\tilde{f}_{j}\\}_{j=1}^{p},\\{\tilde{h}_{uv}\\}_{1\leq
u<v\leq p}\right)$.

###### Proof.

We first prove the result for $p=2$ and then generalize to higher dimensions,
under both sets of constraints on the main effects. Suppose the following
holds for all $x_{1},x_{2}\in[0,1]$:

$\alpha+f_{1}(x_{1})+f_{2}(x_{2})+h_{12}(x_{1},x_{2})=\tilde{\alpha}+\tilde{f_{1}}(x_{1})+\tilde{f_{2}}(x_{2})+\tilde{h}_{12}(x_{1},x_{2}).$
(15)

Proof for origin constraint: We first let $x_{1}=x_{2}=0$. Due to the
constraints on the main effects and interactions, we obtain
$\alpha=\tilde{\alpha}$. Next, we let $x_{2}=0$ for any $x_{1}$. We then
obtain $\alpha+f_{1}(x_{1})=\tilde{\alpha}+\tilde{f_{1}}(x_{1})$ for all
$x_{1}$, which leads us to $f_{1}=\tilde{f_{1}}$. Similarly, by letting
$x_{1}=0$ for any $x_{2}$, we obtain $f_{2}=\tilde{f_{2}}$. Substituting these
results into (15), we find $h_{12}=\tilde{h}_{12}.$

For $p\geq 3$, we fix $1\leq u<v\leq p$. We first let $x_{j}=0$ for all
$j\notin\\{u,v\\}$, which implies

$\alpha+f_{u}(x_{u})+f_{v}(x_{v})+h_{uv}(x_{u},x_{v})=\tilde{\alpha}+\tilde{f}_{u}(x_{u})+\tilde{f}_{v}(x_{v})+\tilde{h}_{uv}(x_{u},x_{v}).$

The result for $p=2$ implies $f_{u}=\tilde{f}_{u}$, $f_{v}=\tilde{f}_{v}$, and
$h_{uv}=\tilde{h}_{uv}$ for any $1\leq u<v\leq p$.

Proof for integral constraint: We first let $x_{2}=0$, which implies
$\alpha+f_{1}(x_{1})+f_{2}(0)=\tilde{\alpha}+\tilde{f}_{1}(x_{1})+\tilde{f_{2}}(0)$
for all $x_{1}\in[0,1]$. Integrating both sides with respect to $x_{1}$ and
using the fact that
$\int_{0}^{1}f_{1}(x_{1})\,dx_{1}=\int_{0}^{1}\tilde{f}_{1}(x_{1})\,dx_{1}=0$,
we obtain $\alpha+f_{2}(0)=\tilde{\alpha}+\tilde{f}_{2}(0)$, which when
substituted into the previous equation, implies that $f_{1}=\tilde{f}_{1}$.
Similarly, letting $x_{1}=0$ leads us to $f_{2}=\tilde{f}_{2}$, which implies
$\alpha=\tilde{\alpha}$ since
$\alpha+f_{2}(0)=\tilde{\alpha}+\tilde{f}_{2}(0)$. Substituting these results
into (15) implies $h_{12}=\tilde{h}_{12}$.

For $p\geq 3$, we fix $1\leq u<v\leq p$. We first let $x_{j}=0$ for all
$j\notin\\{u,v\\}$, which implies

$\rho+f_{u}(x_{u})+f_{v}(x_{v})+h_{uv}(x_{u},x_{v})=\tilde{\rho}+\tilde{f}_{u}(x_{u})+\tilde{f}_{v}(x_{v})+\tilde{h}_{uv}(x_{u},x_{v}),$

where $\rho=\alpha+\sum_{j\notin\\{u,v\\}}f_{j}(0)$ and
$\tilde{\rho}=\tilde{\alpha}+\sum_{j\notin\\{u,v\\}}\tilde{f}_{j}(0)$. Using
the result for $p=2$, this immediately implies that $f_{u}=\tilde{f}_{u}$ for
all $u=1,\ldots,p$, $h_{uv}=\tilde{h}_{uv}$ for all $1\leq u<v\leq p$, and
$\rho=\tilde{\rho}$, which implies $\alpha=\tilde{\alpha}$.

∎

##  Appendix D Other Approaches

We also considered alternative Bayesian monotone spline formulations based on
restricting the sign of the coefficients (Ramsay,, 1988; Shively et al.,,
2009) instead of squaring unconstrained functions. However, squaring
unconstrained functions led to superior estimation performance in simulations,
particularly when the non-negative function being estimated was close to $0$.
To see why, suppose $g(x)=\\{\sum_{j=1}^{m}s_{j}(x)c_{j}\\}^{2}$, where
$\mathbf{s}(x)=(s_{j}(x))_{j=1}^{m}$ denotes $m$ B-spline basis functions and
$\mathbf{c}=(c_{1},\ldots,c_{m})^{\text{T}}$ denotes the vector of basis
coefficients. If we assume $\mathbf{c}\sim N(0,C)$ for some covariance matrix
$C$, then for each $x$,
$g(x)/\left[\mathbf{s}(x)^{\text{T}}\,C\,\mathbf{s}(x)\right]\sim\chi_{1}^{2}$.
Since the $\chi_{1}^{2}$ density function has an infinite spike at $0$, the
squared model provides superior shrinkage of the non-negative function towards
$0$ as compared to the model $g(x)=\mathbf{s}(x)^{\text{T}}\mathbf{c}$, where
we assume the truncated Gaussian prior on $\mathbf{c}\sim
N(0,C)\mathbbm{1}[c_{1}\geq 0,\ldots,c_{m}\geq 0]$.

As an alternative to decomposing $P_{uv}$ and $N_{uv}$ as products of
univariate non-negative functions, we tried modeling $P_{uv}$ and $N_{uv}$ as
squares of linear combinations of tensor products of B-splines (De Boor,,
1978). That is, we let
$P_{uv}(x_{u},x_{v})=\left\\{\displaystyle\sum_{j=1}^{m}\sum_{j^{\prime}=1}^{m}s_{uj}(x_{u})s_{vj^{\prime}}(x_{v})\beta_{uv,jj^{\prime}}\right\\}^{2}$
and
$N_{uv}(x_{u},x_{v})=\left\\{\displaystyle\sum_{j=1}^{m}\sum_{j^{\prime}=1}^{m}s_{uj}(x_{u})s_{vj^{\prime}}(x_{v})\delta_{uv,jj^{\prime}}\right\\}^{2}$.
Assuming $P_{uv}$ and $N_{uv}$ to be products of univariate non-negative
functions is a special case of the above model, obtained by letting
$\beta_{uv,jj^{\prime}}=\theta_{uv,1j}\phi_{uv,1j^{\prime}}$ and
$\delta_{uv,jj^{\prime}}=\theta_{uv,2j}\phi_{uv,2j^{\prime}}.$ This is
equivalent to assuming the matrices
$\mathcal{B}_{uv}=(\beta_{uv,jj^{\prime}})_{1\leq j,j^{\prime}\leq m}$ and
$\mathcal{D}_{uv}=(\delta_{uv,jj^{\prime}})_{1\leq j,j^{\prime}\leq m}$ are
rank $1$. Due to the limited sample sizes and signal-to-noise ratios typical
of environmental epidemiology, we found restricting the rank to $1$ to be much
more effective at identifying synergistic, antagonistic, or null interactions.
Furthermore, the rank $1$ assumption reduces the number of spline coefficients
for the $uv$th interaction from $2m^{2}$ to $4m$, leading to substantial
computational gains.

##  Appendix E Posterior Sampling

Following Section 3.5, we now provide further details for sampling from the
posterior under the proposed SAID approach. Let
$w_{uv}=(\tau_{uv,1},\tau_{uv,2})$. The conditional prior for the spline
coefficients $\mathbf{\Psi}_{uv}$ given the penalty parameter $\kappa_{uv}$
and variance parameters $w_{uv},\nu$ is given by

$\pi(\mathbf{\Psi}_{uv}\mid\kappa_{uv},w_{uv},\nu)\propto\pi_{0}(\mathbf{\Psi}_{uv}\mid
w_{uv},\nu)\exp\left\\{-\kappa_{uv}\widetilde{\mathcal{Q}}(\mathbf{\Psi}_{uv})\right\\},$

where the unpenalized prior $\pi_{0}(\cdot\mid w_{uv},\nu)$ is the product
density given by:

$\displaystyle\pi_{0}(\mathbf{\Psi}_{uv}\mid
w_{uv},\nu)=\displaystyle\prod_{l=1}^{2}\left\\{N(\theta_{uv,l}\mid
0,\nu^{2}\tau_{uv,l}^{2}\Sigma_{0})\,N(\phi_{uv,l}\mid
0,\nu^{2}\tau_{uv,l}^{2}\Sigma_{0})\right\\}.$

Let
$g(\mathbf{\Psi}_{uv},\kappa_{uv},w_{uv},\nu)=\pi_{0}(\mathbf{\Psi}_{uv}\mid
w_{uv},\nu)\,\exp\\{-\kappa_{uv}\widetilde{Q}(\mathbf{\Psi}_{uv})\\}$.
Removing the proportionality, we get:

$\pi(\mathbf{\Psi}_{uv}\mid\kappa_{uv},w_{uv},\nu)=\dfrac{g(\mathbf{\Psi}_{uv},\kappa_{uv},w_{uv},\nu)}{Z(\kappa_{uv},w_{uv},\nu)}=\dfrac{\pi_{0}(\mathbf{\Psi}_{uv}\mid
w_{uv},\nu)\exp\left\\{-\kappa_{uv}\widetilde{\mathcal{Q}}(\mathbf{\Psi}_{uv})\right\\}}{Z(\kappa_{uv},w_{uv},\nu)},$

where the normalizing constant is given by

$Z(\kappa_{uv},w_{uv},\nu)=\int
g(\mathbf{\Psi},\kappa_{uv},w_{uv},\nu)\,d\mathbf{\Psi}.$

With $\pi$ as the prior for $\mathbf{\Psi}_{uv}$, it is difficult to sample
from the posterior of $\mathbf{\Psi}_{uv}$ using standard MCMC algorithms,
since $Z(\kappa_{uv},w_{uv},\nu)$ is an intractable integral. To make MCMC
sampling easier, we follow the variable augmentation technique described in
Rao et al., (2016). Although Rao et al., (2016) develops their algorithm in
the case of intractable likelihoods, we use the same principle to facilitate
joint sampling of
$\mathbf{\Theta}_{uv}=(\mathbf{\Psi}_{uv}^{\text{T}},\kappa_{uv},w_{uv})^{\text{T}}$
and then conditionally sampling $\nu$ given all the parameters
$\\{\mathbf{\Theta}_{uv}\\}_{1\leq u<v\leq p}$, by augmenting new variables
generated using rejection sampling.

We first remark that $g(\mathbf{\Psi}_{uv},\kappa_{uv},w_{uv},\nu)\leq
M\,\pi_{0}(\mathbf{\Psi}_{uv}\mid w_{uv},\nu)$ with the constant $M=1$.
Following Rao et al., (2016), we introduce auxiliary variables
$\mathcal{R}_{uv}=\\{\mathcal{Y}_{uv,1},\ldots,\mathcal{Y}_{uv,|\mathcal{R}_{uv}|}\\}$,
which are considered as the rejected samples before $\mathbf{\Psi}_{uv}$ is
accepted. The rejected proposals $\mathcal{R}_{uv}$ are sampled using the
following rejection sampling algorithm: given $\mathbf{\Theta}_{uv}$ and
$\nu$, we independently sample $\mathcal{Y}_{uv,j}\sim\pi_{0}(\cdot\mid
w_{uv},\nu)$ until the first acceptance; accepting $\mathcal{Y}_{uv,j}$ with
probability
$g(\mathcal{Y}_{uv,j},\kappa_{uv},w_{uv},\nu)/[M\pi_{0}(\mathcal{Y}_{uv,j}\mid
w_{uv},\nu)]=\exp\\{-\kappa_{uv}\widetilde{\mathcal{Q}}(\mathcal{Y}_{uv,j})\\}$.
The conditional joint distribution of $\mathbf{\Psi}_{uv}$ and the rejected
proposals $\mathcal{R}_{uv}$, given $\kappa_{uv},w_{uv}$, and $\nu$ is

$\displaystyle\pi_{1}(\mathbf{\Psi}_{uv},\mathcal{R}_{uv}\mid\kappa_{uv},w_{uv},\nu)$
$\displaystyle=$
$\displaystyle\,g(\mathbf{\Psi}_{uv},\kappa_{uv},w_{uv},\nu)\prod_{j=1}^{|\mathcal{R}_{uv}|}\\{\pi_{0}(\mathcal{Y}_{uv,j}\mid
w_{uv},\nu)-g(\mathcal{Y}_{uv,j},\kappa_{uv},w_{uv},\nu)\\}$ $\displaystyle=$
$\displaystyle\,g(\mathbf{\Psi}_{uv},\kappa_{uv},w_{uv},\nu)\prod_{j=1}^{|\mathcal{R}_{uv}|}\pi_{0}(\mathcal{Y}_{uv,j}\mid
w_{uv},\nu)\left[1-\exp\\{-\kappa_{uv}\widetilde{\mathcal{Q}}(\mathcal{Y}_{uv,j})\\}\right].$

The above joint distribution allows us to compute the prior density of the
$uv$th interaction specific parameters $\mathbf{\Theta}_{uv}$, conditional on
$\nu$ and the rejected proposals $\mathcal{R}_{uv}$:

$\displaystyle
p(\mathbf{\Psi}_{uv},\kappa_{uv},w_{uv}\mid\mathcal{R}_{uv},\nu)$
$\displaystyle\propto\pi_{1}(\mathbf{\Psi}_{uv},\mathcal{R}_{uv}\mid\kappa_{uv},w_{uv},\nu)\,\pi_{\kappa}(\kappa_{uv})\,\pi_{w}(w_{uv}).$
(16)

Here, $\pi_{\kappa}$ is the standard lognormal density; that is, the density
of $\log(W)$ when $W\sim N(0,1)$, and $\pi_{w}(w_{1},w_{2})=C^{+}(w_{1}\mid
0,1)\,C^{+}(w_{2}\mid 0,1)$. From Section 3.5, the conditional likelihood for
$\mathbf{\Theta}_{uv}$ when sampling the $uv$-th interaction is given by
$\mathcal{L}(\mathbf{\Theta}_{uv})=N(\mathbf{\Delta}_{uv}\mid\mathbf{h}_{uv},\sigma^{2}\mathbbm{I}_{n})$.
To sample all the parameters using MCMC, we adopt the following HMC-within-
Gibbs strategy:

  1. 1.

Sample the intercept, main effect, and covariate coefficient parameters as in
Section 3.5.

  2. 2.

For $1\leq u<v\leq p$,

    1. (a)

Construct $\mathbf{\Delta}_{uv}$ after estimating the intercept, main effect,
covariate coefficients, and the rest of the interactions
$\\{h_{u^{\prime}v^{\prime}}:(u^{\prime},v^{\prime})\neq(u,v)\\}$.

    2. (b)

Sample $\mathcal{R}_{uv}$ given $\mathbf{\Theta}_{uv},\nu$: Given
$\mathbf{\Theta}_{uv}$ and $\nu$, we independently sample
$\mathcal{Y}_{uv,j}\sim\pi_{0}(\cdot\mid w_{uv},\nu)$ for $j=1,2,\ldots$ until
a sample is accepted; we accept a sample $\mathcal{Y}_{uv,j}$ with probability
$\exp\\{-\kappa_{uv}\widetilde{\mathcal{Q}}(\mathcal{Y}_{uv,j})\\}$. Once we
accept the first sample, we discard the accepted sample and form the set of
rejected samples
$\mathcal{R}_{uv}=\\{\mathcal{Y}_{uv,1},\ldots,\mathcal{Y}_{uv,|\mathcal{Y}_{uv}|}\\}$.

    3. (c)

Sample $\mathbf{\Theta}_{uv}$ given $\mathcal{R}_{uv}$ and $\nu$: Given
$\mathcal{R}_{uv},\nu,\mathbf{\Delta}_{uv}$ and (16), we form

$\displaystyle\Pi(\mathbf{\Psi}_{uv},\kappa_{uv},w_{uv}\mid\mathbf{\Delta}_{uv},\mathcal{R}_{uv},\nu)\propto\,$
$\displaystyle
N(\mathbf{\Delta}_{uv}\mid\mathbf{h}_{uv},\sigma^{2}\mathbbm{I}_{n})\,p(\mathbf{\Psi}_{uv},\kappa_{uv},w_{uv}\mid\mathcal{R}_{uv},\nu),$

and target sampling from $\Pi$ using HMC.

  3. 3.

Sample $\nu$ given $\\{\mathbf{\Theta}_{uv}\\}_{1\leq u<v\leq p}$ and
$\\{\mathcal{R}_{uv}\\}_{1\leq u<v\leq p}$: Sample $\nu$ from

$\displaystyle p\left(\nu\mid\\{\mathbf{\Theta}_{uv}\\}_{1\leq u<v\leq
p},\\{\mathcal{R}_{uv}\\}_{1\leq u<v\leq p}\right)$
$\displaystyle\propto\displaystyle$
$\displaystyle\,\pi_{\nu}(\nu)\,\prod_{1\leq u<v\leq
p}\left\\{\pi_{0}(\mathbf{\Psi}_{uv}\mid
w_{uv},\nu)\prod_{j=1}^{|\mathcal{R}_{uv}|}\pi_{0}(\mathcal{Y}_{uv,j}\mid
w_{uv},\nu)\right\\}.$

  4. 4.

Sample $\sigma^{2}$ as in Section 3.5.

Since $\nu\sim C^{+}(0,1)$, Step 2 above can be easily accomplished by
introducing an auxiliary parameter $W$ such that $\nu^{2}\mid
W\sim\mbox{IG}(1/2,1/W)$ and $W\sim\mbox{IG}(1/2,1)$ (Makalic and Schmidt,,
2015) and then obtaining conditionally conjugate inverse-gamma updates for
both $\nu^{2}$ and $W$. We now repeat the above steps for a large number of
iterations and discard a burn-in to obtain posterior samples of the unknown
parameters.

##  Appendix F Cutoff in Variable Selection

We provide case 1 and case 2 probabilities for classification using SAID in
Section 4.3, by varying the integral cutoff in $\\{0.005,0.01,0.05,0.10\\}$.

| Cutoff = 0.005 | Cutoff = 0.01 | Cutoff = 0.05 | Cutoff = 0.10
---|---|---|---|---
Synergistic Case 1 Probability | 0.014 | 0.008 | 0.003 | 0.002
Synergistic Case2 Probability | 0.001 | 0.001 | 0.010 | 0.120
Antagonistic Case1 Probability | 0.008 | 0.005 | 0.002 | 0.001
Antagonistic Case2 Probability | 0.010 | 0.001 | 0.010 | 0.200
Null Case1 Probability | 0.001 | 0.001 | 0.010 | 0.160
Null Case2 Probability | 0.024 | 0.014 | 0.005 | 0.001
Table 4: Comparison of different cutoffs when fitting SAID to assess variable
selection accuracy for $p=10$.

##  Appendix G Further tables on simulation results for $p=2$

We provide further tables for the results described in Section 4.2. Tables 5
and 6 provide RMSE values for estimating the dose response surface $H$ and the
interaction $h_{12}$, respectively, under the QR setup. Tables 7 and 8 provide
RMSE values for estimating the dose response surface $H$ and the interaction
$h_{12}$, respectively, under the MIS setup. Details on the QR and MIS setups
are provided in Section 4.2.

Signal and Noise | SAID | BKMR | MixSelect | HierNet | FAMILY | PIE | RAMP
---|---|---|---|---|---|---|---
$\gamma_{0}=1,\sigma_{0}^{2}=0.1$ | 0.04 | 0.03 | 0.03 | 0.17 | 0.43 | 0.21 | 0.03
$\gamma_{0}=1,\sigma_{0}^{2}=0.5$ | 0.10 | 0.07 | 0.06 | 0.17 | 0.44 | 0.24 | 0.08
$\gamma_{0}=2,\sigma_{0}^{2}=0.1$ | 0.05 | 0.04 | 0.04 | 0.17 | 0.57 | 0.62 | 0.03
$\gamma_{0}=2,\sigma_{0}^{2}=0.5$ | 0.10 | 0.08 | 0.06 | 0.18 | 0.56 | 0.62 | 0.07
Table 5: RMSE of competing methods in estimating $H$ for scenario QR. Here $n=500$ and $p=2$. Signal and Noise | SAID | MixSelect | HierNet | FAMILY | PIE | RAMP
---|---|---|---|---|---|---
$\gamma_{0}=1,\sigma_{0}^{2}=0.1$ | 0.08 | 0.08 | 0.33 | 0.31 | 0.05 | 0.04
$\gamma_{0}=1,\sigma_{0}^{2}=0.5$ | 0.12 | 0.11 | 0.29 | 0.31 | 0.10 | 0.17
$\gamma_{0}=2,\sigma_{0}^{2}=0.1$ | 0.12 | 0.08 | 0.37 | 0.62 | 0.05 | 0.04
$\gamma_{0}=2,\sigma_{0}^{2}=0.5$ | 0.17 | 0.09 | 0.35 | 0.62 | 0.11 | 0.10
Table 6: RMSE of competing methods in estimating $h_{12}$ for scenario QR. Here $n=500$ and $p=2$. Signal and Noise | SAID | BKMR | MixSelect | HierNet | PIE | RAMP
---|---|---|---|---|---|---
$\gamma_{0}=1,\sigma_{0}^{2}=0.1$ | 0.06 | 0.05 | 0.10 | 0.19 | 1.44 | 0.08
$\gamma_{0}=1,\sigma_{0}^{2}=0.5$ | 0.12 | 0.12 | 0.12 | 0.20 | 1.10 | 0.13
$\gamma_{0}=2,\sigma_{0}^{2}=0.1$ | 0.01 | 0.02 | 0.04 | 0.05 | 0.52 | 0.03
$\gamma_{0}=2,\sigma_{0}^{2}=0.5$ | 0.08 | 0.08 | 0.14 | 0.19 | 1.59 | 0.11
Table 7: RMSE of competing methods in estimating $H$ for scenario MIS. Here $n=500$ and $p=2$. Signal and Noise | SAID | MixSelect | HierNet | PIE | RAMP
---|---|---|---|---|---
$\gamma_{0}=1,\sigma_{0}^{2}=0.1$ | 0.11 | 0.38 | 0.14 | 0.30 | 0.30
$\gamma_{0}=1,\sigma_{0}^{2}=0.5$ | 0.12 | 0.22 | 0.14 | 0.28 | 0.24
$\gamma_{0}=2,\sigma_{0}^{2}=0.1$ | 0.04 | 0.14 | 0.06 | 0.13 | 0.13
$\gamma_{0}=2,\sigma_{0}^{2}=0.5$ | 0.14 | 0.39 | 0.21 | 0.42 | 0.42
Table 8: RMSE of competing methods in estimating $h_{12}$ for scenario MIS.
Here $n=500$ and $p=2$.

##  Appendix H Creatinine Data Application

We provide plots for rest of the main effects and interaction effects not
presented in Section 5.3. In Figure 6, we plot the estimated main effects of
the metals Antimony, Barium, Cobalt, Lead, Manganese, Strontium, Thallium,
Tin, and Tungsten. In Figure 7, we plot the estimated interaction surface
between Molybdenum and Tin.

|
---|---
|
|
|
|
Figure 6: Plots showing the main effects of dilution-adjusted Antimony,
Barium, Cobalt, Lead, Manganese, Strontium, Thallium, Tin, and Tungsten on log
dilution-adjusted Creatinine. Exposure levels are in log scale. Black line
denotes posterior mean and shaded regions denote pointwise $95\%$ posterior
credible intervals.  Figure 7: Plot showing the interaction effect of
dilution-adjusted Molybdenum and Tin on dilution-adjusted log Creatinine.
Exposure levels are in log scale. Plot shows the pointwise $2.5\%$ posterior
credible surface, the posterior mean, and the $97.5\%$ posterior credible
surface from left to right.