{'fraction_non_alphanumeric': 0.02753368482718219, 'fraction_numerical': 0.012302284710017574, 'mean_word_length': 4.855630126771066, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Flexible Instrumental Variable Models With Bayesian Additive Regression Trees

Methods utilizing instrumental variables have been a fundamental statistical approach to estimation in the presence of unmeasured confounding, usually occurring in non-randomized observational data common to fields such as economics and public health. However, such methods usually make constricting linearity and additivity assumptions that are inapplicable to the complex modeling challenges of today. The growing body of observational data being collected will necessitate flexible regression modeling while also being able to control for confounding using instrumental variables. Therefore, this article presents a nonlinear instrumental variable regression model based on Bayesian regression tree ensembles to estimate such relationships, including interactions, in the presence of confounding. One exciting application of this method is to use genetic variants as instruments, known as Mendelian randomization. Body mass index is one factor that is hypothesized to have a nonlinear relationship with cardiovascular risk factors such as blood pressure while interacting with age. Heterogeneity in patient characteristics such as age could be clinically interesting from a precision medicine perspective where individualized treatment is emphasized. We present our flexible Bayesian instrumental variable regression tree method with an example from the UK Biobank where body mass index is related to blood pressure using genetic variants as the instruments.


Introduction
The presence of unmeasured confounding, particularly in non-randomized observational data, can bias the estimated effect of an exposure on an outcome and make the interpretation of scientific results difficult. Obtaining unbiased results in such a case has traditionally been done through instrumental variable (IV) methods. IV methods incorporate an additional variable, the instrument, to induce exposure variability that is independent of the confounders thereby yielding unbiased estimation of the exposure effect (Stock & Trebbi, 2003). However, this strategy is only valid if the instrument satisfies three properties, usually known as the instrumental variable assumptions. Fields such as economics and public health collect vast quantities of observational data with confounding issues and so researchers in these fields frequently turn to IV analyses to obtain unbiased results.
The three properties central to obtaining unbiased inference are usually known as the instrumental variable assumptions. First, the instruments must be predictive of the exposure in order to induce exposure variability. Second, the instruments must be uncorrelated with the confounders which ensures that this induced variability is independent of the confounders. Finally, there must be no direct effect between the instruments and the outcome, i.e. the instruments only affect the outcome through the exposure. Taken together, these properties imply that associating the instrument-induced exposure to the outcome will be free of confounder influence, yielding an unbiased estimate of the exposure effect that can be interpreted causally. Unfortunately, only the first of these assumptions is testable with the gathered data. As such, it is a notoriously difficult and challenging problem to find valid instruments for a particular research question and it usually requires domain knowledge about the data generating mechanism. For example, a classical question in economics seeks to correctly identify the extent to which education levels correlate with wages, but that relationship is hypothesized to be confounded by the innate ability levels of the subjects (Card, 1999). Not only is ability level unmeasured, but there may not even be consensus as to what constitutes innate ability. Nevertheless, IV methods can be used to assess this research question because ability level only needs to be defined enough to justify the choice of instrument and whether it satisfies the IV assumptions. The geographic proximity to a two or four year college has been used as an instrument in this case, the choice of which is based on domain knowledge within the economics and education literature.
In biomedical research, IV methods have played a large and increasingly important role in the analysis of observational data which is plagued by confounding. Furthermore, the volume and complexity of these datasets are growing over time, necessitating the development of more sophisticated IV methods. See Baiocchi et al. (2014) for an overview of IV methods in biomedical sciences. For example, the use of IV methods has recently shown promise in the area of statistical genetics where a causal genetic relationship can be established by using genetic variants as instruments, specifically single-nucleotide polymorphisms (SNPs), for various traits or phenotypes. This is commonly called Mendelian randomization (MR) because genetic instruments are used as a proxy for the true randomization that is absent in observational data. The IV assumptions in such a case are biologically plausible as long as these SNPs are carefully chosen so as to not be associated with environmental confounders.
Most IV approaches are based on the linear two-stage least-squares (2SLS) methodology which predicts the exposure from the instruments in a first stage regression and then uses the predicted exposure to estimate the unbiased effect of the exposure on the outcome in a second stage regression. Based on the IV assumptions, the variability in the predicted exposure should not be associated with the confounders and so the estimated effect of the exposure will be unbiased. Therefore, the results from this second stage regression will have a causal interpretation assuming the validity of the instruments. However, all 2SLS methods rely on linearity assumptions that can be problematic (Horowitz, 2011). This is particularly true in MR methods that use genetic variants as instruments. Work by Grinberg & Wallace (2021) suggests that linear models may not be sufficient to capture the effect of genetic variants on different traits. Relaxing the linearity assumption could lead to improved prediction in the first stage and improved inference via better power in the second stage. Therefore, methods that relax linearity and additivity assumptions may prove beneficial over ones that do not in the statistical genetics context.
Beyond genetic variants, the traits of interest themselves may have relationships hypothesized to be nonlinear, the estimates of which can provide a more complete picture of the underlying regression relationship. There are a variety of examples in biomedical research (Peter et al., 2015;Scarneciu et al., 2017), including include body mass index with blood pressure, an example that is explored in Section 4 using the UK Biobank dataset. There are additional examples where nonlinear effects have been investigated in fields including, but not limited to: economics (Burke et al., 2015;Botev et al., 2019), marketing (Tuu & Olsen, 2010;Yin et al., 2017), and political science (Kiel, 2000;Lipsitz & Padilla, 2021).
Homogeneity of the exposure effect has been called an implicit assumption of IV methods (Lousdal, 2018), an assumption increasingly regarded as untenable in clinical settings whose focus is on precision medicine and individualized treatment. Individualizing the estimation of the exposure effect by allowing flexible interactions with patient characteristics could be incredibly useful in forming beneficial treatment strategies. Estimating the global effect of a trait without considering age may give an incomplete picture of the results. For example, age is known to have a heterogenous impact on genetic risk for various diseases (Jiang et al., 2021). The National Academy of Medicine lists "increasing evidence generation" as one of five main challenges facing researchers in the area of precision medicine (Dzau & Ginsburg, 2016).
To this end, nonparametric IV regression, which relaxes the linearity and additivity assumptions of 2SLS, has been a popular research area of research for statisticians and econometricians. Examples of such methodologies include Newey & Powell (2003) which are based on series approximations, Burgess et al. (2014) which stratifies the exposure and calculates the causal effect within each strata, and Guo & Small (2016) which adjusts for the first stage errors in the second stage model, something known as control function estimation. In another example, (Chetverikov & Wilhelm, 2017) impose monotonicity which can help stabilize the high variance typically found with series approximations. Stratifying the exposure is an effective strategy but requires the cut points to be chosen by the analyst. This method also utilizes the ratio estimator which can only have a single instrument. Finally, these methods do not allow for general heterogeneity in the exposure effect unless the interactions are specified by the analyst.
Bayesian inference has also been used for IV models as an alternative to methods based on 2SLS or control functions. For examples, see Lopes & Polson (2014) and Wiesenfarth et al. (2014). One benefit is the inherent regularization that is provided by the prior which can help guard against the bias that arises when the instruments are not predictive of the exposure, a violation of the IV assumptions. Another benefit is that both stages are modeled simultaneously with a full probability distribution. This means that the uncertainty inherent in the first stage will be propagated to the second stage during estimation, in contrast to 2SLS where the uncertainty associated with the first stage predictions is not accounted for. Indeed, Bayesian inference allow for the uncertainty for any quantity of interest to be estimated when MCMC is employed to sample from the posterior distribution. One final benefit is that Bayesian inference allows for easy extendibility through hierarchical modeling in order to handle more complicated analysis scenarios.
To allow for both nonlinearity and heterogeneity, Bayesian additive regression trees (Chipman et al., 2010), or BART, priors can be incorporated to the Bayesian framework for IV analysis. BART has shown promise in flexibly estimating general regression relationships without needing any assumptions on the functional form such as higher-order or interaction terms. This article introduces such a heterogenous nonparametric model for use in instrumental variable analyses called npivBART-h, along with simple default settings for the prior so that the method can be used easily and efficiently with minimal learning curve. BART has experienced a usage increase over the past decade along with the development of multiple extensions to handle increasingly complex data. For example, BART has support for many alternative outcome types such as binary, count (Murray, 2021), and survival (Sparapani et al., 2016) including competing risks (Sparapani et al., 2020a) and recurrent events (Sparapani et al., 2020b). Repeated measure outcomes can also be handled through a random effect specification (Tan et al., 2018; that is combined with a BART prior to estimate nonlinear fixed effects. There has also been research toward applying BART in high-dimensional sparse situations (Ročková & van der Pas, 2017;Linero, 2018). Finally, work has been done on applying BART to precision medicine using individidualized treatment rules (Logan et al., 2019). All of these extensions could easily be incorporated into the npivBART-h framework, broadening the scope of its applicability for causal inference.
This article is organized as follows. In Section 2, the npivBART-h model is specified with advice for setting the priors. Section 3 presents a simulation study showing the unbiased inference for these methods in the presence of confounding. Estimation consequences in the presence of weak instruments, an IV assumption violation, is explored through simulation as well. The results from the UK Biobank data relating body mass index to blood pressure with heterogeneity in age and sex using npivBART-h are presented in Section 4. The paper concludes with a discussion in Section 5.

Instrumental Variable Analysis w/ BART
The IV model using Bayesian Additive Regression Trees is described in this section. A brief overview of parametric Bayesian linear IV analysis methods is presented and then our extension to nonlinear relationships using BART is given in Section 2.1. The default regularization priors used in traditional BART are discussed as well as how these priors can be adapted to the IV setting in Section 2.2. Section 2.3 gives a brief treatment of simpler semiparametric models that can estimate interpretable linear effects while Section 2.4 delves into the Dirichlet process mixture model specification of the errors.

Linear Simultaneous Equations Model and npivBART-h
The instrumental variables strategy commonly used in Bayesian IV analyses defines a linear simultaneous equations model as where β is the causal effect of interest assuming the IV assumptions are met. The dependency between equations are adjusted for through the error specification i = ( t i , y i ) ∼ N2(02, Σ). The covariance between the errors controls the degree of confounding. Equation (1) represents the first stage model relating instruments zi with the exposure of interest ti. Equation (2) represents the second stage model with the outcome yi being modeled as a linear function of the exposure and some other measured covariates xi. We assume centered outcomes and so the intercepts are omitted for simplicity. This model is defined in Rossi et al. (2012) and Lopes & Polson (2014). The goal is to relax the linearity assumption in the above model while still retaining unbiased estimation of the causal exposure effect. This can be done with: where i is defined as above. The model that Equations (3) and (4) specify is analogous to the above model, but f2 replaces β as the causal exposure effect of interest. Note that the other covariates xi are incorporated into function f2 which yields a heterogenous causal effect that is not necessarily linear. BART is used to place such a nonlinear prior on f1 and f2, which implies the following specification: Here, g are recursively defined piecewise constant functions otherwise known as regression trees. In this definition, T h and S h represent the probabilistic structure for partitioning the predictor space in function g. M h and L h represent the terminal node values that serve as the output of function g. For simplicity let (T , M) denote the set of all Ht trees approximating f1 and let (S, L) be denoted similarly for f2. In a Bayesian sense, specifying priors on f1 and f2 can be done by specifying priors on (T , M) and (S, L) respectively.

BART and npivBART-h Regularization Priors
All aspects of the BART prior involves some degree of regularization which is what helps guard against overfitting. In particular, Chipman et al. (2010) provide default prior settings that are lightly driven by the data, but perform very well without any cross-validation, although using cross-validation can improve performance. Let us first consider the traditional BART regression model where the target of inference is E[yi|xi] = f (xi) and the regression relationship f is approximated by H regression trees in the ensemble (T , M). For the prior, the trees in this ensemble are taken to be independent. Then, for each (T h , M h ) pair, the prior is broken down into components T h and M h |T h . Partitioning the pair in this way is required because number of terminal nodes in M h is not known until T h is determined.
The prior on the individual regression trees T h involves three components: the probability of a node being an interior node at depth δ, the probability of a splitting rule at an interior node using predictor j, and the probability of the splitting rule at an interior node using cut point c of predictor j. The first is set to decay as the trees grow deeper, favoring simpler trees which helps to avoid overfitting and provides regularization. The second component by default is set to be p −1 for all predictors j = 1, . . . , p. The third component is set by creating an equidistant grid with equal probability for each predictor. This can be extended to binary predictors without loss of generality and categorical / ordinal predictors with a sufficient transformation to dummy indicators. Full details are found in Chipman et al. (2010).
Given the tree structure, the terminal node values M h |T h also need a prior. Let the elements of M h be ν hb for b = 1, . . . , B h where B h is the number of terminal nodes in tree h. A normal prior centered around the sample mean (or 0 if the outcome is centered). The variance of this prior is set so that a certain percentage of the prior captures the range of the outcomes max(yi) − min(yi), reflecting the fact that predictions should not be made outside of this observed range with a high probability. This regularization is quantified by k which controls the prior variability in the the end predictions. For example, k = 2 implies that the range of the data covers the central 95% of the prior while k = 3 implies that this prior coverage is 99.9%.
While these prior settings use the data in the form of the sample mean and range, it should be noted that using the sample mean is regularizing the regression relationship to a flat line which represents the null model in a regression setting. Using the range ensures that the end predictions from our model are calibrated to reasonable values. In the npivBART-h setting, the first stage model can be scaled to the range of the exposure max(ti) − min(ti) and, as in traditional BART, the outcome max(yi) − min(yi). Other components of the BART prior can be modified as context dictates, but that will be unnecessary for the rest of this article.

Semiparametric ivBART Models
Less flexible models are sometimes worth considering because of their added interpretability. For example, suppose that f2(ti, xi) = βti + f22(xi) or alternatively f2(ti, xi) = f21(xi)ti + f22(xi). The former (ivBART-g) estimates a global linear exposure effect after accounting for the nonlinear effect of the measured confounders. The latter (ivBART-h) estimates a linear exposure effect that can vary between observations based on their measured characteristics xi. In particular, the latter estimate could be useful in the biomedical sciences where precision medicine is a major focus of research. Setting the variance of the prior for β in ivBART-g can be done by in the same manner as the previous subsection, but using {max(yi) − min(yi)}/{max(ti) − min(ti)} instead of the range of the outcomes. This modification reflects the fact that β represents a slope term. Method ivBART-g has been discussed in McCulloch et al. (2022).
Finally, npivBART-h has a global counterpart (npivBART-g) that is specified as f2(ti, xi) = f21(ti) + f22(xi) in contrast to the heterogenous version (npivBART-h) introduced in Section 2.1. The priors for npivBART-g can be set in the same way as in the previous subsection, but applied to both f21 and f22.
Estimation for all of these models is trivial using a similar Bayesian backfitting approach outlined in the Supplementary Material for npivBART-h. However, the estimation of ivBART-h can be simplified by using a single regression tree ensemble for both f21(xi) and f22(xi) with a slight modification. Whereas usually a constant mean value is estimated in the terminal nodes M h , here a line with the exposure as the independent variable is estimated. The intercept and slope of this line correspond to f21(xi) and f22(xi). This method is, in principal, equivalent to semiparametric BART models that are particularly useful in a precision medicine context, as done in Zeldow et al. (2019).

Flexible Errors Using Dirichlet Process Mixture Models
As noted in Wiesenfarth et al. (2014) and ignoring the measured confounders xi without loss of generality, 12 t i which is linear in t i . Therefore, the consequence of a bivariate normality assumption is that the model assumes the unmeasured confounders are related to the outcome linearly. This is an unsatisfying assumption when the rest of the model is designed for flexible nonparametric estimation. As such, this assumption can be relaxed using a Dirichlet process mixture (DPM) model. The DPM model can be thought of as an infinite mixture of normals that can model any probability distribution and so the induced linearity assumption is relaxed with such a specification. An additional benefit of this specification is that it is robust to outliers because it can model an error distribution with fatter tails than the usual normal errors.
This flexible error specification can be written as with Σi ∼ F and F ∼ DP(F0, α) which defines a DPM model. F0 represents the base distribution for Σi and α is the sparsity parameter. In this case, a two-dimensional inverse-Wishart distribution is used as F0. As detailed in the Supplementary Materials, sampling can be done using the algorithm of Escobar & West (1995).

Simulation
The simulation study aims to answer two separate issues concerning nonparametric instrumental variable regression using BART: (1) the overall effectiveness of using BART ensembles for an IV analysis and (2) the effect of using a Bayesian analysis for IV in the presence of weak instruments. Results from each simulation are presented in this section.
To judge the effectiveness of the proposed method, data under eight varying simulation settings are generated 1,000 times repeatedly. The sample size is set as n = 1,000 and the variance of the errors is set as V ( t i ) = V ( y i ) = 1. This implies that Cov( t i , y i ) = ρ, the correlation. We set ρ = 0.0 and ρ = 0.7 to account for the situation where there are no unmeasured confounders and the situation where there are unmeasured confounders respectively. The regression function for the first stage is f1(zi) = sin(πz1iz2i) + 2z 2 3i + z4i + 0.5z5i is the traditional Friedman function (Friedman, 1991), commonly used for benchmarking nonlinear prediction problems. The results from transformation f1 are scaled by their variance so that there is an equal amount of signal and noise in stage 1. Vector zi was simulated according to the allele frequencies and correlation structure of the SNPs found in gene FTO, a gene associated with BMI and other cardiovascular related diseases. SNPs from this gene are used as the instruments for the causal analysis of UK Biobank data in Section 4. 20 SNPs are simulated for zi. Ten additional independent covariates were simulated as xi ∼ Unif(−1, 1). The functional form of the true second stage regression relationship is varied according to Table 1.
Four versions of the method are considered, each corresponding to one of the true regression relationships above. These versions estimate: a global linear exposure effect (ivBART-g), a global nonlinear exposure effect (npivBART-g), a heterogenous linear exposure effect (ivBART-h), and a heterogenous nonlinear exposure effect (npivBART-h). The first three are the models discussed in Section 2.3. Additionally, two-stage least-squares (2SLS) and simple BART (fitting the second stage model without instruments) are also considered as comparators.
For each model, f2(ti, xi) is estimated for ten test points in the covariate space. These points are the five values of t on an equidistant grid between −2.5 and 2.5. These five points are then used in conjunction with x1 = −0.5 and x1 = 0.5 (representing the x1 < 0 and x1 ≥ 0 cases respectively). The partial dependence function is used to margin out the xj for j = 2, . . . , p. This yields ten estimates of f2 which are compared to the true value of f2 at each of these points. Then, the bias can be computed at each grid point while the RMSE can be computed across the five grid points for each of the positive and negative x1. The bias results for t ∈ {0, 1.25, 2.50} are given in Figure 1 and Figure 2 which give the results in the situation with ρ = 0.0 and ρ = 0.7. Table 2 gives the results for the RMSE across the grid points in every scenario. Note that the negative gridpoints {−1.25, −2.50} show similar results, though the direction of the bias is usually flipped.
In particular, it is obvious that traditional BART, although it has as much flexibility as npivBART-h, cannot estimate f2 without bias when ρ = 0.7 as expected. The bias from BART is especially large for grid point t = 2.5 in the sparse tail of the data. The model that is expected to correctly capture the regression relationship performs better than its less flexible counterparts. For example, ivBART-g and ivBART-h perform poorly when the relationship is nonlinear while ivBART-g and npivBART-g perform poorly when the relationship is heterogenous. 2SLS performs poorly with either nonlinearity or heterogeneity.
There does appear to be a slight benefit for choosing a linear model when the underlying truth is linear. Intuitively, a regression coefficient is more stable than estimating an ensemble of regression trees. More flexible nonlinear models still provide unbiased estimation in the linear case, but have an increased estimation variability. Table 2 presents the RMSE results with their standard errors and confirms that the linear methods have slightly reduced error than the nonlinear methods when the truth is linear, though both are unbiased.
An additional simulation was used to assess how weak instruments, one violation of the IV assumptions, may affect estimation. For simplicity, only the global linear variation of f2(ti, xi) is considered as the true regression   Traditional BART is biased when t = 1.25 or t = 2.5 implying that the confounding bias is worse in the tails where the data is sparse. The mean (sd) RMSE for estimates of f 2 across the grid of t values and 1,000 simulations which mimic the bias results. However, it is of note that the linear methods perform slightly better than the nonlinear methods when the underlying truth is linear. This is because estimating linear coefficients is more stable than tree-based nonparametric estimators and this reduced variability yields reduced error (bias being equal). relationship. Likewise, only ivBART-g is considered so that the simulation criteria is the resulting estimated value of β. To mimic the weak instrument scenario, the true f1(zi) was scaled by a constant C ∈ {0.0, 0.1, 0.5, 1.0}. A larger sample size of n = 5,000 was considered in addition to the sample size of n = 1,000 used in the first simulation.
The scalar k represents the number of standard deviations covering the range of the data on the prior for the end predictions. This is described in Section 2.2 as well as Chipman et al. (2010). For example, the default k = 2 implies that central 95% of the prior covers the range of the data while k = 3 implies that the central 99.9% of the prior covers the range of the data. Hence, increasing k results in stronger regularization of the predictions (via stronger regularization of β) to 0. The results in Table 3 confirm that weak instruments will result in biased estimation. The BART-based methods slightly improve estimation of β compared to 2SLS likely due to the regularization, but both are biased when C is small. As expected, a large enough sample size could overcome this bias as indicated by the cases with small C and n = 5,000. Table 3: The mean (sd) estimate of β over 1,000 replications from the weak instrument scenarios when β = 1. Estimation appears to be biased in the weak instruments cases as expected, i.e. when C = 0.0 or C = 0.1. The value k controls the amount of regularization through the amount of prior probability placed around 0. A higher value of k means that the resulting posterior distribution will be closer to 0. Additionally, increasing the sample size results in less biased estimation as expected. n = 1,000 n = 5,000

IV Analysis Using UKBB Data
The large scale UK Biobank (UKBB) data which has approximately 500,000 volunteers enrolled for study between the ages of 40 and 69 serves as a demonstration of the potential utility of this method. The centralized UK National Health Service allows for all health related events including disease diagnoses, drug prescriptions, and deaths to be recorded. Beginning in 2017, researchers were able to access the genetic information recorded with the database including genetic variants such as SNPs. As such, this database can be used to relate body mass index (BMI) as an exposure with systolic blood pressure (sBP) using single nucleotide polymorphisms (SNPs) as instruments.
BMI and blood pressure (systolic and diastolic) are traits known to have a strong relationship (Landi et al., 2018) that is thought to be nonlinear. Additionally, age and sex are two covariates that could be associated Table 4: Summary statistics for sBP and dBP by BMI, sex, and age. These summary results suggest an effect of all three on sBP. There does not appear to be a strong effect of either age or sex when looking at dBP as an outcome, however. As such, only sBP will be considered in the subsequent heterogenous analysis. with either exposure or outcome. It would be interesting to explore whether the causal relationships differ based on either of these patient characteristics. These questions can be answered using the npivBART-h methodology introduced in this article.
To start, the summary statistics for systolic blood pressure (sBP) and diastolic blood pressure (dBP) outcomes by different groupings of BMI, age and sex are displayed in Table 4. The summary statistics show a fairly large difference in the sBP distribution between the groups of BMI and age with a smaller difference between sexes. However, the summary statistics of dBP do not imply much of distributional difference in either age or sex. Therefore, dBP is not considered as an outcome for the rest of this section. Of course, these summary statistics do not speak to any potential interaction between the predictors, but that can be addressed when looking at the npivBART-h results.
Gene FTO on chromosome 16 is known to be associated with obesity related disorders and contains the first genetic variant to be associated with BMI. Discovered in 2007, subsequent studies have discovered additional genetic variants from FTO that are also strongly associated with BMI and this effect appears to be consistent among different ethnic populations. SNPs from the cis-region of this gene, defined as 50kbp on each side of the gene, of FTO were extracted as candidate instruments for this analysis. Of these extracted genes, pZ = 8 SNPs were used as instruments because they have been identified as risk variants for BMI in a genome-wide replication study (Tan et al., 2014). Such genetic variants should not be associated with environmental confounding factors, motivating their use as instruments and is commonly called Mendelian randomization in the statistical genetics literature.
In total, 473,303 observations were used from the UKBB dataset. The Dirichlet Process Mixture (DPM), defined in Equation (5), was used as a flexible error distribution. Three models are used for comparison in this section. The first is an IV regression that is nonlinear in BMI using npivBART-h, the second is a nonlinear regression without IVs using traditional BART, and the third is an IV regression that is linear in BMI using ivBART-h. For all methods, there were 5,000 burn-in samples and 5,000 additional samples drawn from the posterior distribution across three chains leading to a total of 15,000 posterior samples. Each chain utilized different starting values to promote proper mixing. Thinning the MCMC sample by a factor of k to reduce the autocorrelation, while recommended in the past, is now recognized as not being beneficial because reducing the total number of posterior draws by a factor of k is more detrimental than autocorrelation in most cases (Link & Eaton, 2012). As such, no thinning is performed for these analyses.
The resulting curves given in Figure 3 were plotted for a grid of fifteen points of BMI, four ages, and both sexes. The nonlinear IV method (left) utilizes a nonlinear second stage for the IV regression. First, it appears that there may be a three-way interaction between BMI, age, and sex. For the male observations in particular, the effect of age is much larger for moderate values of BMI compared to large values of BMI where the age effect is reduced. This two-way interaction does not appear in the female group leading to the three-way interaction. When comparing the nonlinear IV method to its non-IV counterpart (middle), there appears to be a difference in the estimated curves, particularly in the tail values of BMI. The non-IV method tends to estimate a steeper increasing trend overall. This seems reasonable because the non-IV method is confounded by environmental factors, making the effect of BMI appear stronger than it truly is. Finally, the linear IV method (right) is considered as comparison to the nonlinear methods. This plot suggests a much stronger three-way interaction than its nonlinear comparators. The rigidity of linear estimates are apparent here because the predictions for higher BMI values are overwhelmed by the observations with moderate BMI values where the increasing effect is strongest. In cases such as this, nonlinear regression methods will paint a more complete and nuanced picture of the true relationship.
While a difference between the IV and non-IV methods appears visually, it can be confirmed by looking at the estimated error correlation. The relationship is confounded by unmeasured variables when this correlation is non-zero. However, because the DPM is used as a prior for the error distribution, each observation has its own 2 × 2 variance-covariance matrix estimated at each draw of the MCMC chain. Let the correlation term from this variance-covariance matrix be denoted as ρ id for observation i at MCMC iteration d. Let averages of these Figure 4: The distribution of posterior means across the data (left) and the distribution of data means across the posterior (right) for the error correlation ρ. The three chains are indicated with shading. The left plot shows little heterogeneity within the error correlation among the observations. As such, averaging over the sample as displayed in the right yields a posterior distribution that indicates ρ > 0. estimates be defined asρ ·d = (1/n) n i=1 ρ id andρi· = (1/D) D d=1 ρ id where D = 5,000 in this case. Note that ρi· are computed within each chain for visualization. Histograms of these distributions are given in Figure 4 with each of the three chains distinguished using the shading. The left panel shows little correlation heterogeneity among the observations within each chain, suggesting that a global estimate of this correlation is sufficient. In fact, most iterations of the MCMC chain resulted in less than 3 clusters being detected among the observations. Because of this, averaging over the observations to yield the posterior distribution of a single correlation term will be valid. The right panel of the figure, which shows this posterior distribution, indicates that the error correlation is positive and non-zero, though relatively small. This is in line with the relatively small visual difference seen between the IV and non-IV curves of Figure 3.

Discussion
This article introduces npivBART-h, a highly flexible instrumental variable framework using Bayesian additive regression trees capable of estimating the individualized nonlinear relationship between an exposure and outcome in the presence of unmeasured confounding. Simulation results confirm that this is a feasible method for estimating heterogenous nonlinear regression relationships with a causal interpretation. With the increasingly large volume of complex observational data being gathered today, statistical methods to flexibly estimate the regression relationship while also accounting for unmeasured confounding will necessary. These nonlinear estimates can give researchers a more thorough set of results. Indeed, there are a variety of applications where applying flexible methods such as npivBART-h could prove scientifically interesting. This article presents one such application in the area of statistical genetics and Mendelian randomization where genetic variants are used as instruments to compute a causal effect of traits on an outcome. As done in Section 4, body mass index can be causally related to blood pressure, a relationship that is commonly thought to be nonlinear.
Additionally, the method allows for individualized heterogenous estimation which could prove useful for tailoring treatments to the characteristics of the patient, usually known as precision medicine. This is an exciting new frontier in the medical field and could lead to individualized treatment strategies that improve health outcomes for society at large. Developing analysis tools such as npivBART-h to execute this type of research is an opportunity for statisticians to meaningfully contribute to this potential paradigm shift.
As mentioned in Section 1, the Bayesian probability framework has many benefits. First, it allows for easy extendibility through hierarchical models so that more complex data can be analyzed. One simple extension to our method would be to allow more complicated outcomes such as binary, count, or time to event either for the exposure or outcome. These types of outcomes have already been considered for traditional BART and so this extension would be a relatively straightforward way to broaden the method's applicability. Additionally, sparsity based methods such as Ročková & van der Pas (2017) and Linero (2018) should be considered to exclude potential weak instruments in the first stage which can bias the results as shown through simulation. While it may boost the predictive power of the model, including additional instruments increases the likelihood of IV assumption violations. It is for this reason that a parsimonious set of SNPs is to be preferred, something that can be achieved using sparsity methods.
One interesting application in statistical genetics would be to apply this methodology across an entire genome, relating the expression levels of a gene to an outcome in a possibly nonlinear fashion. In this case, the expression levels of each gene use the SNPs within their cis-region as instruments. Doing this across the entire genome is known as a transcriptome-wide association study (TWAS) and is usually based on linear 2SLS methods. However, there is precedent for using machine learning prediction models for this type of analysis (Okoro et al., 2021). For an overview of TWAS and additional statistical considerations, see Xue et al. (2020) or Li & Ritchie (2021).
An R software package entitled hnpivbart is available to fit any of these models. This implementation of the BART algorithm is based off of the BART3 package  which utilizes C++ code for fast and efficient operation of the method.