{'fraction_non_alphanumeric': 0.025849704813048265, 'fraction_numerical': 0.006820986624862413, 'mean_word_length': 5.041003425347572, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonparametric Estimation of Population Average Dose-Response Curves using Entropy Balancing Weights for Continuous Exposures

Weighted estimators are commonly used for estimating exposure effects in observational settings to establish causal relations. These estimators have a long history of development when the exposure of interest is binary and where the weights are typically functions of an estimated propensity score. Recent developments in optimization-based estimators for constructing weights in binary exposure settings, such as those based on entropy balancing, have shown more promise in estimating treatment effects than those methods that focus on the direct estimation of the propensity score using likelihood-based methods. This paper explores recent developments of entropy balancing methods to continuous exposure settings and the estimation of population dose-response curves using nonparametric estimation combined with entropy balancing weights, focusing on factors that would be important to applied researchers in medical or health services research. The methods developed here are applied to data from a study assessing the effect of non-randomized components of an evidence-based substance use treatment program on emotional and substance use clinical outcomes.


Introduction
Assessing causal relationships between continuous treatments and clinical outcomes is often of interest in substance use research. For example, we may seek to understand how substance use outcomes are affected by the dose of a substance use treatment program. Evaluations of such programs can be analytically challenging because researchers are often unable to control the factors (observed and unobserved) that would allow them to fully understand the relationship between the treatment program and the clinical outcomes. The gold standard for making causal claims is a randomized experiment in which control over important variables would be attained through randomization. Use of observational studies in health care domains, such as in our motivating example assessing the effect of non-randomized components of an evidence-based substance use treatment program on clinical outcomes, offer more difficult settings with which to assess causal relationships. Strong assumptions must be made to defend claims that would be considered even close to causal in these settings. Statistical methods from the field of causal inference provide much needed tools that enable researchers to come closer to addressing the questions they would hope to address using a clinical trial, assuming such strong assumptions can be asserted.
The field of causal inference has an extended history for addressing the challenges inherent in observational studies, often focused on methods involving the propensity score (Rosenbaum & Rubin 1983) and settings where the exposure of interest is binary (see, Rosenbaum (2010), Imbens & Rubin (2015) for examples and extended history). More recent extensions of these methods to settings where the exposure of interest is not binary have arisen, and are often focused on the estimation of population dose-response curves (Robins et al. 2000, Imbens 2000, Hirano & Imbens 2004, Imai & van Dyk 2004 that describe the relationship between a continuous treatment and an outcome. Many of these extensions in observational data settings estimate the dose-response curve conditional on the "generalized propensity score" (GPS) through stratification or inverse probability weighting (IPW). Unfortunately, these methods appear to be under-utilized in the medical and health care literature, potentially because there is not enough guidance for the research community to readily adopt these methods and/or because their use requires asserting strong assumptions that may be unlikely to hold in practice. An assumption of the correct specification of the propensity score is difficult enough to assert in the binary setting; the assumption of the correct distribution of the GPS needed in estimating the dose-response curve in the continuous treatment setting can be even more difficult.
There has been some work to improve the utility of methods for continuous treatments or exposures by focusing on flexible modeling of the GPS and outcome models (Zhu et al. 2015, Kennedy et al. 2017. Unfortunately, flexible models can sometimes be insufficient to fully reduce bias in observational studies when models are misspecified (e.g., the GPS requires modeling the entire distribution of the exposure given observed covariates and misspecification can occur if flexible modeling is used to only estimate the mean relationship between exposure and covariates). It may be that there is more promise, as has been shown in the binary settings, when using optimization based weighting methods that seek to directly optimize "covariate balance" (Hainmueller 2012, Imai & Ratkovic 2014, Zubizarreta 2015, Vegetabile et al. 2020. These methods remove the requirement of modeling the propensity score or GPS directly. One such method, entropy balancing, has recently been demonstrated to have properties similar to doubly robust estimators (Zhao & Percival 2017), where if either the outcome model or assignment model is specified correctly, then consistent estimation is possible. These advances in optimal balancing have begun to be applied to continuous exposure settings (Yiu & Su 2018, Fong et al. 2018, Kallus & Santacatterina 2019 and recent work has extended the entropy balancing methodology to this setting (Tübbicke 2020). These methods hold great promise for the estimation of more accurate doseresponse curves in health care and medical research, and a richer understanding of dose-response curves has the potential to save lives and improve patient outcomes.
In this paper, we examine the use of entropy balancing for continuous treatments and illustrate how it can be used with nonparametric regression methods for estimating population dose-response curves. Our work is distinct from concurrent research on this topic (Tübbicke 2020) in the following ways: (1) we consider nonparametric methods for the estimation of the dose-response curve that better reflect real world relationships; (2) we provide a careful discussion of the influence that higher order moments of the exposure and covariates have on the entropy balancing algorithm and subsequent estimation of the dose-response curve; (3) we provide evidence on the performance of bootstrap-based confidence intervals for point-wise inference of the dose-response curve; and (4) we include a user-friendly R package for deriving entropy balancing weights. Throughout the paper, we emphasize considerations that must be made in practice in order to provide researchers with needed guidance on implementing these methods. Our motivating case study examines the causal effects of the number of treatment components received for an evidence based substance use treatment program for adolescents on both substance use and emotional well-being outcomes. The use of nonparametric estimation and entropy balancing together removes the need for the researcher to choose a correct model specification of the GPS and a correct form of the outcome model, removing two key sources of potential error in analyses.
The organization of the paper is as follows: Section 2 provides the development of the methodology, including an overview of causal inference in the potential outcomes framework for continuous exposures, the development of entropy balancing for continuous exposures, careful consideration of key issues for applied analysis (i.e., the assessment of covariate balance and choosing the number of moments to balance), and details on inference, the local linear regression strategy, and bootstrapped standard error estimators. Section 3 provides a simulation that attempts to stress the method in reasonable ways (e.g., non-linearity in the exposure/covariate relationships and non-linear outcome models). Section 4 applies the methods to data collected from over 2200 adolescents substance users receiving the Adolescent Community Reinforcement Approach (A-CRA) (Godley et al. 2016) treatment program with the goal of assessing the causal effect of the number of behavioral treatment components received as part of their A-CRA program, on longitudinal substance use and emotional well-being outcomes. Section 5 provides a concluding discussion and remarks on the use of these methods for health services researchers.

Methods
In this section, we provide an overview of the methods for nonparametric estimation of population average dose-response curves using entropy balancing. We begin with an overview of causal inference in the continuous exposure setting that introduces our notation throughout. We then provide an extension of entropy balancing for continuous exposure settings and demonstrate how these weights can be used with local regression estimation techniques that can incorporate the entropy balancing weights and nonparametrically estimate dose-response curves. As stated in the introduction, we have developed these methods into an R package (project development can be found at https://github.com/bvegetabile/entbal).

Notation and Causal Inference Framework
We utilize a generalization of the potential outcomes framework of Neyman et al. (1923) and Rubin (1974), specifically focusing on notation from the extension to general exposure settings outlined in Hirano & Imbens (2004). We define an index set i ∈ {1, . . . , N} of units for that we are interested in assessing some "dose-response" relationship, i.e., a functional relationship between an exposure and an outcome within some population. We define the exposure variable for each unit as A i taking values from the set A, e.g., A ∈ + , and define a random vector X i ∈ X as the set of variables that occur prior to the application of the exposure A i and are responsible for setting its value; typically X i = (X i1 , X i2 , . . . , X ij , . . . , X id ) will be a d-dimensional vector and embedded in d . We also define the observed outcome of interest as a variable Y i and we posit that for each unit there exists a potential outcome Under this framework, an ideal measurement process would observe an infinite number of copies of each unit, apply the exposure at exactly the same instance in space and time across these copies and measure each outcome Y i (a); this is clearly not possible in most practical cases and in the case where there are only two exposures, this has been referred to as the Fundamental Problem of Causal Inference (Holland 1986). Therefore, we leverage information across units and estimate the functional f (a) = E[Y (a)] ≡ X Y (a)p(x)∂x, where the expectation is with respect to a population of individuals defined by the random vector X. The functional f (·) defines the population dose-response curve, i.e., the expected response to the application of A = a in a population defined by the vector X. Using this functional, we could choose to define causal effects as being contrasts of f (a), e.g., There are numerous estimation strategies for the estimation of f (·) within the Neyman/Rubin framework, e.g., regression estimation (Hirano & Imbens 2004) or through the use of subclassification on the GPS (Imai & van Dyk 2004). In this paper, we utilize weighted estimation similar to the marginal structural approach (Robins et al. 2000). The core of the procedure is the result that if we define where p(·) denotes a density function, then we can obtain consistent estimation for the average response, i.e., it follows that Note that the denominator of Equation 1, i.e., e(a, x) ≡ p(a|x), is referred to as the GPS in the literature and the weights are referred to as stabilized inverse probability of treatment weights (IPTW).
The results in Equation 2 require a few strong assumptions (see Imai & van Dyk (2004) and Hirano & Imbens (2004) for overviews) for consistent estimation and a causal interpretation, specifically: 1) an extension of the stable unit treatment value assumption (SUTVA) for continuous exposure settings, and 2) ignorability assumptions. Assumption 1 states that the distributions of the potential outcomes for a unit are assumed independent of the potential exposure status of another unit, given observed covariates. This assumption implies that there is no interference between units and that we can use a consistency argument to assert that the observed outcome under treatment A = a is equal to the corresponding potential outcome, i.e., Y i = Y i (A i ). The ignorability assumptions contain two requirements: 1) unconfoundedness, and 2) positivity. Unconfoundedness requires that the potential outcomes are conditionally independent of treatment given the covariates, e.g., a weak unconfoundedness assumption defined in Hirano & Imbens (2004) states that Y (a) ⊥ ⊥ A | X for all a ∈ A. Ignorability assumptions also require positivity, i.e., 0 < p(A ∈ T |X) for all X ∈ X and measurable sets T ⊂ A with positive measure (Imai & van Dyk 2004). The implication of the ignorability assumptions are that all covariates responsible for assigning exposure levels have been identified and any missing or unobserved X * variable is a potential source of estimation bias. They also imply that for any given unit with X = x there needs to be sufficiently high probability of receiving any treatment level under consideration. These are both strong assumptions, e.g., unconfoundedness is untestable (Imbens & Rubin 2015), that require thought and defense in any analysis. Additionally, positivity may not hold if it is observed that there is a strong relationship between any covariate and the exposure. For this paper, we assume that these assumptions hold.
Using these assumptions, we will use weighted estimation methods to estimate the population dose-response curve and develop a method for constructing optimal weights.

Entropy Balancing for Continuous Exposures
Entropy balancing (Hainmueller 2012) is a method of creating optimal weights for causal estimation. Within the causal literature, the method was originally developed for the estimation of causal effects in binary exposure settings, though related methods have a longer history in survey calibration (Haberman 1984, Deville & Särndal 1992, Deville et al. 1993).
Here we demonstrate an extension of entropy balancing to handle continuous exposures. This has been proposed in Tübbicke (2020), but we adopt a slightly different notation to provide guidance on the influence that higher moments of the exposure and covariates have in the entropy balancing algorithm and the subsequent estimation of dose-response curves. If we define µ p j = E[X p j ], the mean of an element X j from vector X, and µ q A = E[A q ] for integers p and q and assume that we have the correct weights as defined in (1), then the method relies on the observation that in the joint distribution of A and X, the weighted covariance between X p j and A q for any pair (p, q) should be zero, while all other weighted marginal relationships should remain the same. This is demonstrated for the covariance between X and A here, This observation will allow us to select nonparametric weights, i.e., those that do not model p(a) or p(a|x) directly, that enforce the constraints that we would expect to observe if the weights had the correct ratio-ofdensities interpretation. A similar constraint was utilized in Fong et al. (2018) for their nonparametric covariate balancing propensity score (npCBPS) for generalized treatment regimes.
To select weights, we will minimize the Kullback-Leibler divergence between two distributions of weights: one that utilizes the weights as defined above and another based on empirical weights, i.e., q i ≡ N −1 , subject to constraints on weighted covariance between A and X and enforcing constraints on the marginal distributions. This minimization has the effect of finding a set of weights w that are "closest" in the sense of KL-divergence to a set of weights that maximize the sample information. We will also impose additional constraints that ensure the weights are positive and sum to one. This is formalized below: This can be optimized using the method of Lagrange multipliers, where the objective becomes The solution to this objective is of a similar form to that found in Hainmueller (2012), but here the weights are of the form, .
Plugging these weights into Equation (4) we find the dual form that is no longer a function of the weights, This new dual objective can now be optimized using an efficient convex optimization routine and has a reduced parameter space. Our implementations have found the limited memory Broyden-Fletcher-Goldfarb-Shanno algorithm with bounding constraints (L-BFGS-B) implemented in R to be computationally efficient. The above optimization is in some cases numerically unstable, i.e., computational overflows can occur when exponentiating the sum of the multiplications of λ, γ, φ and powers of each observation's constraint variables (e.g., X p ij − µ p j ), and therefore we further constrain the optimization to be bounded such that parameters lie in the range λ, γ, φ ∈ [l, u]. In practice, these are set to large values, e.g., [l, u] = [−100, 100] and can be further constrained if there are numerical instabilities. Note that in practice when these numerical issues occur, there is often an issue in satisfying the optimization constraints. The bounding constraints are used in our implementation to ensure that the algorithm at least returns a set of weights to the user that can be evaluated. This behavior though implies that "perfect" balance with respect to the constraints is not possible in this data, which may be observed by remaining imbalances when performing covariate balance assessments, or if the λ, γ, φ are equal to the constraint values. When these parameters must be constrained to lie in a small range in practice, it should be assessed how sensitive the results are to this bounding choice.
For more details on the original formulation of this problem in binary settings, see Hainmueller (2012) and Appendix A for an overview.

Choosing the number of moments to balance
The above algorithm requires specifying the number of moments, i.e., the parameters (p, q), that must be balanced. As suggested by Equation (3), the definition of the weights implies that the marginal distributions of A and X should remain unchanged under weighting and that only the covariances should be zero. Clearly though, as more constraints are placed on the optimization procedure it is more challenging to find a solution set and typically the effect of additional moment constraints is a reduction in the "effective" number of usable observations for analysis (discussed further in Section 2.2.2). To keep the marginals correct for A and X, we have observed that in our simulations balancing two or three moments when variables are continuous and one moment when they are binary (ordinal variables are treated as continuous and multi-valued variables are converted to binary variables through dummy coding) has provided adequate performance. The higher moment, i.e., p, q = 3, may be helpful in cases where the treatment or covariates are skewed, or if there are non-linearities in the relationship between A and elements of the vector X. Similarly, we do not remove covariances between X p j and A for p > 3, as typically such strong relationships would imply that the positivity assumption is violated, or the effects of higher moments are weak enough that they can be captured with lower moments.

Assessing weighted balance and the effect of weighting on sample information
An important step in any causal analysis is assessing "covariate balance" and the reduction in sample information due to the weighting scheme employed (captured by the effective sample size in analyses that utilize the weights). Covariate balance in the binary setting is typically assessed by comparing both unweighted and weighted standardized mean differences or other test statistics (e.g., Kolmogorov-Smirnov statistics) to understand how well the weights have improved the overall comparability of the two groups. Visual inspection is also often performed comparing the weighted empirical cumulative distribution function (CDF) to the CDF of the target population for inference for the covariates. The idea behind these balance assessments is that when the weights are applied to each group (e.g., exposure versus control) the weighted distribution of the pre-exposure confounders will be similar to the target population over which inference is being estimated (e.g., for the average treatment effect in the population (ATE) the target is the overall population characteristics; for the average treatment effect among those that are treated, the target is the characteristics of the treatment group).
In the continuous exposure setting, we will perform similar diagnostics from the binary setting, but include more comparisons to ensure we have adequately assessed covariate balance. Note that these balance assessments are not specific to entropy balancing weights and would apply for any type of weights that may be used in a continuous exposure setting, but throughout we focus on how they apply in our setting. A comprehensive assessment of covariate balance should focus on the following summaries between the weighted and unweighted distributions: • Conditional mean balance: linearly regressing A on each X j and assessing the model coefficient to provide the magnitude of the relationship between the two variables • Correlation balance: assessing the Pearson correlation between A and each X j • Marginal balance: assessing marginal summary statistics of the variables in the data set, i.e., each X j and A, specifically calculating the following summaries: The conditional mean balance is a similar assessment as in the binary setting, i.e., we are leveraging the fact that if our weights have the correct ratio-of-density interpretations then . This observation combined with our constraint that the covariance must be zero imply that the model coefficient in a weighted regression using our entropy balancing weights should also be zero if the constraints can be satisfied. The conditional mean balance is different than the binary setting where we do not perform a regression and simply compare the weighted means in each group, i.e., comparing E[w X|A = a] and E[w X|A = a ] to see if they are the same.
The next assessments, correlation balance and marginal balance, investigate if the weights have achieved zero correlation between the exposure and each covariate while maintaining the original marginal distributions over which inference is required. One concern would be that if the weighted exposure distribution is different than the unweighted distribution, then it implies that the weights have changed the distribution of the exposure over which inference is possible. To assess correlation balance we leverage the fact that the primary constraint in the optimization was that the weighted covariance between each X p j and A is zero. To fairly compare the performance across covariates on a standardized scale we instead look at the Pearson correlation coefficient between X p j and A in both the weighted and unweighted distributions. If the optimization is successful, then these values should be very small (up to the tolerance of the optimization routine), though it has been suggested elsewhere that when using the GPS, values less than 0.1 should be acceptable in practice (Zhu et al. 2015). The optimization also specifies that certain moments of the marginal distributions of the exposure and covariates are unchanged in the weighted distribution. To validate this, we compare the unweighted and weighted mean and variance of the exposure variable and covariates to assess that they are unchanged by the weighting. Similarly, the KS-statistic provides an additional assessment of the distortion of the marginal distribution and we compare the KS-statistics between the weighted and unweighted distributions. As the algorithm does not balance all moments (typically we set the max (p, q) to 3), we do not expect the KS-statistic to be exactly zero, but they should also not be large. In our experience with entropy balancing for continuous exposures, KS-statistics less than 0.1 are common and typically imply small differences in the tails of the distributions (in cases where the mean and variance have been matched). If any of these values (correlation balance or marginal balance) appear large, then an optimal solution may not be possible and/or the data may not be appropriate for estimating a population dose-response curve.
Finally, we perform an assessment of the effective sample size as employing weights in an analysis typically reduces the sample information relative to an equal weighting; specifically greater variability in the weights corresponds to a greater loss of information. One reason for only focusing on balancing two or three moments is that by adding more moments we place more constraints on the optimization procedure. This often requires increasing the variability in weights and the effect is that this can make estimates of the dose-response relationship more variable. The reduction in sample information is typically measured using the "effective sample size (ESS)" (Kish 1965), a measure from the survey literature related to the concept of a "design effect". The effective sample and in a binary setting would typically be calculated within each group, but in the continuous setting it is calculated using all weights. Within an applied analysis, the effective sample size can be thought of as the number of "usable" units that meet the constraints placed on the sample through weighting, i.e., if a unit's weight is effectively zero relative to all other unit's weights then that observation does not contribute to the analysis. Additionally, checking the loss of information is useful for assessing the impact of weighting and determining if the available data support the analysis goals, i.e., small effective sample sizes between 0 and 100 imply low power and high variability in point estimates (i.e., low power will inflate the standard errors of the estimator).

Estimation of Population Dose-Response Curves
The weights obtained in the previous section can be utilized in the same fashion as survey weights in order to estimate the population average dose-response curve; therefore we focus on flexible estimation techniques that can incorporate weights. An overview of the estimation procedure in this paper is as follows: we utilize local linear regression (Cleveland & Devlin 1988, Friedman et al. 2001) that is further weighted by entropy balancing weights to understand the population average behavior in a neighborhood around an exposure A = a 0 of interest. We can then extend this to estimate the dose-response curve in the high-density region of the distribution of A. We stress that our focus on flexible estimation of the balancing weights and the dose-response imply that we will only be able to estimate within the high density region of the exposure distribution and extrapolation past this high density region is not recommended.
To estimate locally the method requires defining a target point of interest A = a 0 and selecting a kernel function κ(a 0 , ·) that defines the "closeness" between the target point a 0 and any other observed point A i = a. Within this local region the model assumes a linear model where, The coefficients β 0 (a 0 ) and β 1 (a 0 ) are therefore a function of the target point, i.e., they are locally adaptive. Estimation then follows by minimizing a loss function within the neighborhood defined by the kernel function and further weighted by our entropy balancing weights, e.g., if Y is continuous and choosing to minimize squared error loss, The above requires specifying a kernel function κ(a 0 , ·). A common example for local linear regression is the tricube function, κ tc (a 0 , ·) = 1 − |u| 3 3 , where |u| ≤ 1 and each u i is a normalized distance function, i.e., u i = d(a 0 , A i )/ max j d(a 0 , A j ). The distance function d(·, ·) typically contains a nuisance parameter α that controls the fraction of points contained in the neighborhood of points used in estimating the outcome (Cleveland & Devlin 1988). Other examples of kernels include the squared exponential or cosine similarity, but are not considered here (Scholkopf & Smola 2001). It is possible to use cross-validation to select the parameter α, or another model selection procedure, and in our application we use two-fold cross validation as each fold requires a large number of observations to minimize the sensitivity of the out-of-sample estimates across folds to the weight distribution. If necessary, this procedure can be repeated multiple times to assess the distribution of the resulting α parameter.
To implement these methods, we use the loess package in R that implements the tricube kernel and efficient estimation. For an extensive overview of local regression methods in the context of prediction models see Chapter 6 in Friedman et al. (2001).

Inference using Bootstrapped Confidence Intervals
Section 2.3 provides point estimates of the average treatment effect in the population at specific points A = a 0 , but does not provide adequate information for quantifying uncertainty. Therefore to provide uncertainty estimates, we use a bootstrap procedure (Efron & Tibshirani 1994) focused on bootstrapping the entire estimation procedure. Specifically for b = 1, . . . , B, we sample with replacement from the observations to create a new data set D b = {Y b , A b , X b } and using this new sample to re-estimate the entropy balancing weights and perform the entire dose-response estimation procedure on the outcomes Y b to obtain a new estimatef b (A i ) for all i . That is, we perform cross-validation to select a new α b , and refit the local linear model to obtain a new estimate of the treatment effect at that point. This process is repeated for the set of B bootstrap samples and the standard error of the bootstrap distribution of estimates at the point A = a 0 , denoted se(f b (a 0 )), is calculated. Approximate confidence intervals can then be constructed, i.e., We evaluate properties of this interval estimator in Section 3.4.

Data Generation
To demonstrate the method we perform a simulation study that attempts to stress the method in reasonable ways, i.e., non-linear relationship with A in the outcome model and a strictly positive treatment variable. We generate 1000 simulations where for each replication we generate a sample size of N = 1000.
We generate data with three covariates, and relate these to treatment in the following way, where µ A is a non-centrality parameter that creates a somewhat realistic treatment distribution generated from a mixture of covariate types, i.e., high concentration of treatment at zero with a longer right tail of higher valued treatments, where an example of such an exposure may be the cumulative number of procedures that a client is exposed to during a substance abuse treatment program similar to our motivating example.
The outcome is generated as follows, This outcome model introduces non-linearity in A that would be difficult to capture with traditional modeling and includes an interaction between treatment and covariates that can cause a strong bias when only modeling the marginal relationship between A and Y and there is confounding by X. An additional simulation is included in the appendix that considers the case where there is no effect of treatment.
The left panel of Figure 1 demonstrates the marginal distribution of A and the right panel demonstrates the observed relationships between A and Y with the true population marginal dose-response curve overlain and a linear smoother through the data. The figure demonstrates that there is confounding for this population doseresponse curve, where if A and X did not have strong relationships, then the curves would be similar. Additionally, the trend in the unweighted distribution is in the opposite direction, i.e., in the observed relationship there appears to be a positive association between exposure and outcome when the true marginal relation (averaged over the population) is a negative quadratic association.

Estimation Methods Compared
To fairly evaluate the method, we compare entropy balancing against other methods of estimating weights; these are listed in Table 1. We primarily focus on comparing against two "naïve" type analyses and an alternative using the generalized Covariate Balancing Propensity Score (CBPS) of Fong et al. (2018) and implemented in the CBPS R package; specifically their fully nonparametric version. For both entropy balancing and CBPS, we perform the estimation focusing on a range of moments balanced.
To evaluate the performance of the dose-response curve, we propose estimation using the suggested local linear regression method, as well as a parametric "global" Weight Estimation Method 1) Unweighted estimation; i.e., naïve analysis 2) Entropy balancing (1); as described above, one moment -One moment of X covariance balanced (i.e., p = 1) with A (q = 1) -One marginal moment balanced (i.e., p, q = 1) 3) Entropy balancing (2); two moments.

Performance Balancing Covariates
The first result compares correlation balance, specifically the correlation between each covariate and the exposure variable. This is illustrated in Table 2 where the correlation balance statistics for the entropy balancing method are all less than 1e −3 (both the average across replications and the maximum value attained). We also see that the method provides good performance in retaining the marginal distributions under weighting. The KS-statistics across replications are all less than 0.1, which is often considered a good guideline for adequate covariate balance in binary settings. Table 2 also demonstrates that using misspecified linear models can provide unpredictable behavior and erratic covariate imbalance measures, where the average absolute weighted correlation is higher than what was observed in the sample. Additionally, there are large discrepancies between the weighted covariate distributions and the observed distributions. Finally, the CBPS methods perform similar to the entropy balancing method. The CBPS methods obtain higher weighted correlations on average than entropy balancing, but the values would be very acceptable in practice.

Conditional Balance
Correlation Balance Marginal Balance: KS statistics only

Performance in Estimating Dose-Response Curves
While assessing covariate balance is an important step in a causal inference analysis, the weights are often nuisance parameters in the estimation of the effect of interest, i.e., in this application the population dose-response curve. Figures 2 and 3 illustrate visually the performance of the methods in estimating the doseresponse curves, i.e., the sampling distribution across the simulation replications. In these figures, it is clear that the unweighted estimates provide curves that would incorrectly indicate a positive association despite the true marginal relationship being a negative quadratic association. Methods that focus on balancing a higher number of moments (i.e., 2 or 3) appear to capture the true relationship well in this simulation and the average of the curves across the simulations within the high-density region match the true marginal relationship. Additionally, the distribution across the curves demonstrates that there is not a high degree of variability in these estimates despite the effective sample size being reduced to a small fraction of the original number of observations on average. For methods that focus on only balancing one moment, we see that these weights reduce the bias in estimating the population dose-response curve, but that they do not capture the true relationship for high exposure values, i.e., A > 20.   a quantification of the differences. Estimation using local linear regression and balancing on more than one moment provides better performance than balancing only one moment; this finding is observed for both the proposed entropy balancing method and the CBPS. In this example, the entropy balancing method with higher moments generally performs better than the CBPS method when using local linear regression, but both methods are more similar when the true outcome model is known. One surprise is that despite the poor covariate balance, the linear model performed better than CBPS when only one moment is matched in this simulation setting. All methods reduce the MSE relative to an unweighted analysis because they are controlling for confounding variables through weighting.

Coverage of Bootstrap Confidence Intervals
The results of the previous two sections demonstrate that the entropy balancing procedure is both able to provide adequate covariate balance and in some instances better performance in providing point estimates for the population dose-response curve in the high density region of the exposure. A point estimate though is often only part of the inferential story and a rigorous analysis requires a quantification of the uncertainty about that point estimate. To that end, we construct a simulation study that examines the performance of the bootstrap procedure defined in Section 2.4.
In this simulation we use the data-generation procedure outlined in Section 3 and again construct 1000 replications, but for each estimation procedure also provide bootstrap confidence intervals using B = 100 replications. The results from this study are provided in The results demonstrate that the standard errors from the bootstrap estimation procedure are slightly conservative when compared with the sampling distribution of the estimates (lower right panel). This results in confidence intervals that are often too wide and therefore provide greater than 95% coverage as demonstrated in the right panel in the high density region of the exposure. Coverage suffers towards the tails where there may be some bias in estimation due to sparsity of observations in this region. Future work should examine the properties of this procedure for estimating uncertainty in weighted population dose-response curves.

Overview of the case study data and analysis
We provide a demonstration of these methods using data collected from four cohorts of the Assertive Adolescent and Family Treatment Initiative (AAFTI) that had, as part of its initiative, a goal of assessing the effect of the Adolescent Community Reinforcement Approach (A-CRA) (Godley et al. 2016(Godley et al. , 2011 on longitudinal substance use and behavioral outcomes among adolescent substance users. A-CRA is a community-based reinforcement behavioral treatment program for alcohol and substance use disorders administered over a three-month period that recommends 10 individual focused sessions between a clinician and a client, two sessions between the clinician and at least one caregiver, and an additional two sessions that include all parties. However, clinicians are able to suggest as many sessions as necessary within an individual treatment program. Within each session the therapist performs a set of procedures from a therapy guide (outlined in Godley et al. (2016)), such as the "Problem Solving Skills Training" designed to improve the clients ability to solve problems outside of therapy, the "Happiness Scale" used to understand levels of happiness in various aspects of the clients life, and "Increasing Prosocial Recreation" that focuses on developing behaviors that increase interactions with individuals that are less likely to lead to drug use, among many other procedures. While the therapy guide is a suggested course of treatment, each trained clinician administering A-CRA is able to use additional sessions if necessary and has latitude in how they choose the many procedures at their disposal during treatment sessions. Individuals within AAFTI were then also exposed to Assertive Continuing Care, that had A-CRA embedded within its sessions, for an additional three months; providing a total of sixth months of A-CRA exposure. Among those adolescents engaging in treatment, investigating the relationship between the number of sessions an individual receives and emotional and substance use outcomes is of interest. In the AAFTI, clinicians engaged in 4 to 79 sessions, where 4 sessions was a criteria for engagement with the substance use treatment. Individuals that were identified as "not engaging" with treatment were removed from the analysis.
Throughout the implementation of A-CRA, covariate information on the client is collected and additional targeted information related to psychological risk factors and substance use are collected through the Global Appraisal of Individuals Needs (GAIN) (Dennis et al. 2003), a survey composed of measures that attempt to capture domains such as substance use, physical health, mental and emotional health, living situation, and legal situation. Within AAFTI, both covariate information and the GAIN were measured at baseline, 3 months, 6 months, and 12 months, providing longitudinal assessment over a one-year period. In this analysis we focus on assessing the marginal relationship between the number of treatment sessions delivered to the adolescent client on two outcome variables of interest at the 6-month follow-up point: 1) the emotional problems scale (EPS) of the GAIN, and 2) the binary indicator of whether the individual is "in recovery" or not. The variable "in recovery" is defined as being one (true) if the client has been housed in the community, has been abstinent from any substance use within the past month, and reports no past month substance use problems. Figure 5 demonstrates that there is an increase of those in recovery at six months when compared to baseline.

Number of Sessions
In addition to the primary exposure (i.e., number of treatment sessions) and outcome variables, we consider baseline covariate information such as age, gender, race, and a subset of the GAIN scales that have been identified as being necessary to control for in our analysis. Figure 5 contains visualizations of the initial distributions of key variables (exposure and outcomes) to highlight similarities with our simulation study. One feature of the exposure distribution is that the 99th quantile is approximately 46, so we do not expect to trust inference from our method much higher than this value, as we would be extrapolating to a part of the exposure distribution with low representation. We estimate the dose-response curve in the range [4,50] and remove individuals not in this range when performing the analysis. Within the original data set, there is missingness within some of the covariates. In this analysis, we analyse a single imputed data set from a larger imputation procedure of these data used to test the relative effectiveness of A-CRA versus another evidence based treatment program for adolescent substance users (Griffin et al. in review). As we are focused on highlighting and demonstrating the utility of the method, we leave discussions of the full implementation of these methods with missing data to future work. Table 4 contains a detailed list of the covariate information and baseline GAIN measures that were included in this analysis.
There are two complications related to our primary exposure variable that complicate an applied causal analysis in our case study: 1) the exposure variable is not decided at entry, that is, there are factors related to the number of treatment sessions that cannot be observed longitudinally throughout exposure to A-CRA, such as client and clinician interactions, and 2) there can be multiple versions of the same treatment level as we do not know the combination of the underlying procedures used (e.g., two clients each exposed to 10 sessions could have an underlying different set of procedures that they were exposed to). These issues introduce complications in interpreting our analysis as causal that we discuss further Section 4.3. Regardless, the entropy balancing weights should provide a robust associational relationship that marginalizes over, i.e., controls for, the observed covariates and will yield useful insights for clinicians about the relationship between total number of A-CRA treatment sessions and our outcomes. Table 4 organizes important features of the data set, such as summaries of the marginal relationships and correlations of covariates with the primary exposure variable (number of treatment sessions), both without and with weighting using the entropy balancing weights (matched on three-moments). While there are only mild correlations among the covariates and the exposure variable prior to weighting, many of these relationships have regression coefficients that are statistically significant. The table also demonstrates that the entropy balancing weights provide adequate covariate balance as many of the correlations are small and close to zero and there are no longer relationships present when regressing the exposure variable on each covariate in the weighted data. Another important result is that the weighted distribution of each covariate is relatively unchanged under this weighting scheme (due to the marginal constraints of the algorithm) as illustrated through the similar weighted and unweighted mean and standard deviations, as well as the KS-statistic comparing the weighted and unweighted distributions. Finally, there is only a small reduction in effective sample size in performing this analysis, which is not surprising given the mild correlations in the unweighted data. Figure 6 presents the results from performing weighted local linear regressions using the entropy balancing weights. The left panels show the estimated curve and the bootstrapped 95% confidence interval overlain onto the original data and the center panels focus in on the relationships. The right panels provide a closer look at the curves and compares against other methods of estimating the outcome variables. Specifically, we demonstrate the results using a naïve (unweighted) local linear regression, a naïve (unweighted) regression, i.e., Y ∼ A, a regression controlling for baseline variables of interest (EPS or "in recovery" status at baseline), and a weighted regression using the entropy balancing weights and a simple model, i.e., Y ∼ A. Comparing across these results, we see that our method provides results that are somewhat distinct from the other methods of estimation and that the results for the in recovery variable and the EPS are similar (but trending in opposite directions). Combined with the simulation results of the last section, this suggests that the method is able to more accurately estimate the true marginal relationship and that the use of nonparametric estimation frees the  Figure 6: Estimated relationships between the number of treatment sessions and the emotional problem scale (EPS) at the 6-month follow-up (top) and "in recovery" status at the 6-month follow-up (bottom). The right panel of each row compares estimates for the proposed method (EB + LOESS) to 1) a naïve (unweighted) local linear regression, 2) a naïve (unweighted) regression modeled such that exposure is the only covariate, i.e., Y ∼ A, 3) a simple linear regression controlling for the baseline value of the outcome measure, i.e., Y ∼ A + X, and 4) a weighted linear regression only controlling for the exposure variable, i.e., Y ∼ A.  Table 4: Summary statistics of covariate information, both without (left section) and with (right section) weighting using the entropy balancing weights. The last column presents the KS-statistics comparing the weighted and unweighted covariate distribution for each covariate. Each section is organized as follows: mean and standard deviation of the covariate (marginal balance), correlation with the exposure variable (correlation balance), and the regression coefficient and corresponding t-statistic and p-value obtained by regressing exposure on the covariate indicated (conditional mean balance). The last row of the table represents the effective number of observations without and with weighting. The table illutrates that the entropy balancing removes the observed correlation among the covariates and the exposure and that they are now conditionally independent given the weights.

Relationship between Number of Sessions and In Recovery Status at 6−Months
analyst from having to determine the correct functional relationship.
These results also provide a few interesting insights into the relationship between A-CRA and our outcomes. The first result is that there is a positive approximately quadratic relationship between the exposure variable and EPS at 6-months; those with the lowest and highest number of treatment sessions appear to have "worse" outcomes on EPS. Similarly, there is a negative approximately quadratic relationship between the number of treatment sessions and in recovery status; those exposed to fewer sessions do similar to those with many sessions, and both perform worse on this key substance use outcome than those toward the center of the exposure distribution. Taken causally, these results would be confusing because it appears that the number of treatment sessions causes worse outcomes on the EPS and substance use outcomes, when a monotonic relationship may be expected if the number of sessions is playing a role in these outcomes. As was discussed earlier, the number of sessions likely violates some of the key assumptions for a causal interpretation and a more reasonable interpretation is that there is an interaction between the clinician and the client throughout the A-CRA exposure that explains these conflicting results. That is, it could be that individuals receiving fewer sessions are not truly engaging in treatment and this does not provide enough "dose" in the form of number of treatment sessions to the client. Subsequently, they perform worse on the substance use outcome of being in recovery. Similarly, those that have a larger number of sessions are more difficult to treat, and therefore we see worse outcomes at these levels. Compared with those two groups, those in the middle of the exposure distribution are more "typical" cases and we see these individuals having the best outcomes. This essentially suggests a form of unobserved confounding and highlights the need for well-defined treatments and critical thinking about an ideal experimental set-up that one would want to perform to address the scientific question. The lack of a rigorous causal interpretation does not diminish that these are interesting results that provide useful insights to clinicians, such as the fact that more information should be recorded and considered at intake on the individuals and throughout their substance use program to better understand and act upon the underlying differences among patients at different exposure levels.

Discussion
This manuscript outlined a nonparametric estimation procedure using local linear regression that incorporates entropy balancing weights for performing inference in observational settings. The method was demonstrated through a comprehensive simulation study and through an applied example that assessed the relationship between the number of treatment sessions received during substance use treatment on substance use and emotional wellbeing outcomes among adolescent substance users. The simulation study demonstrated that the methodology can provide comparable (or in some cases superior) performance to a competing method (CBPS) and that both out perform naïve analyses using linear models. Additionally, the number of moments that are balanced in the analysis is critical. When only the first moment was balanced, the methods did similar to naïve analyses; including higher moments in the balancing equations further reduces the bias in the analysis for both entropy balancing and the CBPS. The application demonstrated that the method can be useful in health care settings and provide insight to clinicians, while freeing them from modeling decisions in the estimation phase (though it does not free them from critically thinking about assumptions for causal inference and underlying data generation).
Causality is often the underlying intention (albeit in some cases a subconscious intention) when collecting data in an attempt to organize phenomenon in the world around us and requires a fine control on many conditions in the world. Establishing "true" causality is often unattainable though (see the Fundamental Problem of Causal Inference) and we must make assumptions that enable us to find comparisons of outcomes in the world under counterfactual (i.e., alternative) exposures. The method developed here is no different and requires many strong assumptions that can be implausible and in some cases untestable. The method supposes that a suitable estimand in an observational study is the average dose-response effect within a population, an analog from clinical trials, that provides insight on the broader behavior of an outcome as it relates to the exposure. Researchers should be cautioned though that the positivity assumption is necessary to defend this estimand and it requires that for each value X = x there is a positive probability of receiving any exposure level in the range of exposures of interest (something guaranteed in a completely randomized experiment) in order to enable estimation. However, in many observational settings, this can be an unrealistic assumption. Any form of moderate preferable selection of exposures that can be ascribed to covariates, both observed and unobserved, makes satisfying this assumption less plausible. When there are strong relationships in the data, then estimation of the conditional or local average exposure effect should be considered by thinking through the suitability of average effects in some neighborhood of X = x. One diagnostic for assessing the positivity assumption in observational settings is assessing the effective sample size. When this measure tends toward small numbers, there is not enough sample to meet the constraints within the balancing algorithm and implies that the population dose-response curve estimand should not be addressed with the data.
Even when a positivity assumption appears plausible, there are additional hurdles to making causal claims that a researcher must consider such as temporal ordering, the structure of the relationships between covariates and the exposure, unconfoundedness, and model selection. Our methodology attempts to relieve the researcher of one of these (i.e., model selection) by using local linear models combined with entropy balancing, but the researcher must defend and document assertions that the remaining assumptions hold. In particular, in mental health and health care settings, both the structure of the relationships among variables and the unconfoundedness assumption must be made with care. In our example, it would appear that temporal ordering and positivity are plausible, but there appears to be some unobserved variables in the relationship between the number of sessions and the outcomes due to the longitudinal nature of the decision process in "assigning" treatment. These deficiencies do not imply that the application of these methods do not inform science, but they should come with the appropriate caveats that this work mainly informs hypothesis generation and recommendations for future work. A related topic within the field of causal inference are sensitivity analyses that attempt to isolate how sensitive a result is to assumptions made within the design pipeline. Future work should address and develop sensitivity analyses in observational studies with continuous exposures in nonparametric modeling settings.
The methods developed here provide a unique flexibility for estimation in observational settings that incorporates weights to remove biases. While the assumptions required to make a causal claim are numerous, the combination of nonparametric estimation and entropy balancing weights can clearly remove many biases in observational studies and should become a common tool for applied researchers.

A Binary Exposure Entropy Balancing Algorithm -Hainmueller (2012)
Consider an outcome variable Y , a binary exposure variable A, a multivariate vector of covariates X, and a distribution represented by p(x) over which inference is required with specified moments µ p j = E X [X p j ] ≡ X x p j p(x)∂x. The goal in the binary setting is estimating weights w i , w i ≥ 0, i w i = 1 such that the estimator in Equation (6), i.e., differences in potential outcomes, has consistent estimation for Entropy balancing selects the weights such that the weighted empirical moments within each group A = a match those from the population of interest, i.e., To do this it minimizes the KL-divergence between the distribution of the weights that satisfies the desired distributional characteristics (represented by constraints on the weighted moments) and a "base" distribution of weights q (typically defined as q i ≡ N −1 representing the empirical weights) within each group: Hainmueller (2012) chose to optimize this object using the method of Lagrange multipliers, The objective above provides a closed form solution for the weights that can dramatically reduce the dimensionality of the optimization problem and facilitate finding an efficient solution using off-the-shelf optimizers. See Hainmueller (2012) for a full treatment of these methods in the case of binary treatments for estimating the average effect on the exposed.

B Supplemental Simulation -No Effect of Exposure
In this supplemental material, we reconsider the simulation of Section 3, but modify the outcome variable so that there is no direct effect of the exposure variable on the outcome variable Y . To accomplish this, in this simulation the data generation procedure is augmented such that Y = X 1 + X 2 1 + X 2 + X 2 2 + X 1 X 2 + X 3 + , ∼ N (0, 1) and all other data-generating parameters remain the same. This implies that the marginal effect is E[Y (a)] = 3.55 for all a ∈ A and thus A has no population-level effects. Within the simulation we again model using the local regression method, but here compare against a global linear model where the assumption is that Y (a) ≈ β 0 + β 1 A to assess how the weights perform without requiring a selection of the correct form of the outcome model.
A visualization of the induced relationship between exposure and outcome is shown in Figure 7.  Figure 7: Visualization of the distribution of treatment variable (left) and relationship between treatment and outcome (right). In the left figure, the high-density region of A is highlighted. In the right figure, the true marginal relationship is shown in red, a simple unweighted linear smoother is shown in blue, and an unweighted linear estimation is provided in green.

B.1.1 Covariate Balance
The covariate balance results are unchanged as the relationships among A and X in the data-generation are the same and the seed used for the random number generator is the same.

B.1.2 Performance in Estimating Dose-Response Curves
The performance here is similar to the results of Section 3. Both entropy balancing and the CBPS perform well when a high enough number of moments is balanced and all other methods perform less than satisfactory. The bootstrap confidence intervals again remain conservative, but in this case it is observed that at certain points of the exposure distribution there is under-coverage. Future work should further investigate confidence intervals for these methods.