A Boosting Algorithm for Estimating Generalized Propensity Scores with Continuous Treatments

2.1 Definition and assumptions

Let Y denote the response of interest, T be the treatment level and X be a p-dimensional vector of baseline covariates. The observed data can be represented as (Yi,Ti,Xi), i=1,…,n, a random sample from (Y,T,X). In addition to the observed quantities, we further define Yi(t) as the potential outcome if subject i were assigned to treatment-level t. Here, T is a random variable and t is a specific level of T. The dose–response function we are interested in estimating is μ(t)=E[Yi(t)], and we assume Yi(t) is well defined for t∈τ, where τ is a continuous domain.

Similar to the binary case, the ignorability assumption is as follows:

where f(t|⋅) refers to the conditional density. That is, the treatment assignment is conditionally independent of the potential outcomes given the covariates. In other words, we assume there are no unmeasured covariates that may jointly influence the treatment assignment and potential outcomes.

Denote the generalized propensity score as r(t,X)≡fT|X(t|X), which is the conditional density of observing the treatment-level t given the covariates [6]. The ignorability assumption implies

That is, to adjust for confounding, it is sufficient to condition on the generalized propensity scores instead of conditioning on the vector of covariates, which might be high dimensional.

2.2 Estimation based on marginal structural models

Under the ignorability assumption, we focus on the marginal structural model approach to estimate the dose–response function proposed by Robins [9] and Robins et al. [13]. The method works by building a marginal structural model for the potential outcomes. For example, we may assume a linear model:

Model (1) is marginal because it is defined for the expected value of potential outcomes without conditioning on any covariates (which is different from regression models). Based on the observed data (Yi,Ti,Xi), i=1,…,n, the parameters in eq. (1) can be consistently estimated by IPW. For the ith subject, the weight is

There are two important issues related to this approach: (i) the estimation of the inverse probability weights; (ii) the functional form of the outcome model in eq. (1). The first issue is the main topic of this article and will be explored in the next section. Here we briefly discuss the second issue. The consistency result of MSM estimators relies on the correct specification of the outcome model. However, the true form of E[Y(t)] is unknown in reality and a flexible model is always preferred. In the data application, we assume a regression spline function [14] for the outcome model:

where u+p=(u+)p is a truncated power function and u+=max(0,u). p is the order of the polynomial function and K is the total number of inner knots. The inner knots are either distributed evenly on τ or defined as the equally spaced sample quantiles of Ti, i=1,…,n. That is,

where T(i) is the ith quantile of T1, T2, …, Tn.

To determine the regression spline function, we need to find the optimal p and K. The traditional model selection criteria, such as AIC and BIC, are based on a simple random sample. These criteria can be extended to determine the form of the marginal structural model for a non-randomized sample [15, 16]. Under the assumption that Y is normally distributed, the weighted AIC can be defined as

where l is the total number of parameters. In eq. (3), l=K+p+1. Notice that in this stage, we treat wi′s as fixed. Similarly, we define the weighted BIC as

We will illustrate the specification of the outcome model through the data application in Section 5.3. In the next section, we will focus on the first issue and propose a boosting algorithm for estimating the generalized propensity scores.

3.1 Modeling the generalized propensity scores

In the MSM approach, the estimation of wi′s in eq. (2) is essential. For simplicity, we assume T follows a normal distribution so that r(Ti) can be easily estimated by normal density. To be noticed, if the normality assumption does not hold for T, we can always employ a nonparametric method, such as Kernel density estimation, to estimate r(Ti). To estimate r(Ti,Xi), a traditional way is to assume

Then, the estimation of r(Ti,Xi) follows two steps [13]:

1. Run a multiple regression of Ti on Xi, i=1,…,n and get Tˆi and σˆ;

2. Calculate the residuals εˆi=Ti−Tˆi; r(Ti,Xi) can be approximated by

   (7)rˆ(Ti,Xi)≈f(εˆi)≈12πσˆexp−εˆi22σˆ2.

where m(X) is the mean function of T given X. We advocate a machine learning algorithm, boosting, to estimate m(X). The advantage of boosting is that it is a nonparametric algorithm that can automatically pick important covariates, nonlinear terms and interaction terms among covariates [12]. It fits an additive model and each component (base learner) is a regression tree. Mathematically, it can be written as:

where M is the total number of trees, Km is the number of terminal nodes for the mth tree, Rmj is the indicator of rectangular region in the feature space spanned by X and cmj is the predicted constant in region Rmj. Km and Rmj are determined by optimizing some nonparametric information criterion, such as Entropy, misclassification rate or Gini Index. cmj is simply the average value of Ti in the training data that fall in the region Rmj. Details about how to construct a tree classifier can be found in Breiman et al. [17].

In boosting, M is an important tuning parameter. It determines the trade-off between bias and variance of the causal estimator. In inverse weighted methods, if subject i receives a weight wi, it means the subject will be replicated wi−1 times; that is, there will be wi−1 replications in the weighted pseudo-sample. In the weighted sample, if the propensity scores are correctly estimated, the treatment assignment and the covariates are supposed to be unconfounded under the ignorability assumption [13]. Consequently, the causal effect can be estimated as in a simple randomized study without confounding. Therefore, a reasonable criterion is to stop the algorithm at the number of trees such that the treatment assignment and the covariates are independent (unconfounded) in the weighted sample. Next, we propose stopping criteria for boosting based on this idea.

3.2 Algorithm

We propose four different criteria (summarized in Table 1) for how to measure the degree of independence/correlation between the treatment and each covariate. These criteria are named as “Pearson/polyserial,” “Spearman,” “Kendall” and “distance.” Pearson/polyserial, Spearman and Kendall correlations are commonly used correlation matrices; distance correlation [18, 19] is the most recently proposed and is gaining popularity due to its nice property: it can be defined for two variables of arbitrary dimensions and arbitrary types. Next, we will briefly describe these four correlations.

We denote one of the covariates as Xj, for j=1,2,…,J, where J is the total number of covariates. If both T and Xj are normally distributed, the Pearson correlation coefficient will be zero given that T and Xj are independent. When Xj is categorical, the Pearson correlation coefficient could be biased [20]. Instead, we should use the polyserial correlation coefficient [21], which essentially assumes that the categorical variable is obtained by classifying an underlying continuous variable Xj′. The unknown parameters of Xj′ can be estimated by maximum likelihood. Then, the polyserial correlation is calculated as the Pearson correlation between T and Xj′. Spearman and Kendall correlation coefficients are rank-based correlations that can be applied to both continuous and categorical variables. If T and Xj are independent, we would expect the Spearman and Kendall correlation coefficients to be close to zero. A more flexible measurement of correlation/independence is distance correlation. The distance correlation takes values between zero and one and it equals zero if and only if T and Xj are independent, regardless of the type of Xj.

In the binary treatment case, to check whether the propensity scores adequately balance the covariates, we calculate the standardized difference in the weighted mean between the treatment group and the control group. If balance is achieved, the difference should be small. In the continuous case, we propose a general algorithm that uses bootstrapping to calculate the weighted correlation coefficient. The procedure requires the following steps for each value of M (number of trees).

1. Calculate rˆ(Ti,Xi) using boosting with M trees. Then, calculate

   wi=rˆ(Ti)rˆ(Ti,Xi)fori=1,…,n.

2. Sample n observations from the original dataset with replacement. Each data point is sampled with the inverse probability weight obtained from the first step. Calculate the corresponding coefficient between T and Xj on the weighted sample and denote it as dji;

3. Repeat Step (2) k times and get dj1,dj2,…,djk. Calculate the average correlation coefficient between T and Xj, denoted as dˉj;

4. Perform a Fisher transformation on dˉj, i.e.,

   (10)zj=12ln1+dˉj1−dˉj.

5. Average the absolute value of zj over all the covariates and get the average absolute correlation coefficient (AACC).

For each value of M=1,2,…,20,000, calculate AACC and find the optimal number of trees that lead to the smallest AACC value. The R code for calculating AACC is displayed in Appendix A. An alternative suggestion is to replace Step 5 by calculating the maximum value of the absolute correlation coefficient (MACC) over all the covariates and find the optimal number of trees that lead to the smallest MACC value. After the value of M is determined, the generalized propensity score is estimated by eqs (9) and (7). The Fisher transformation in Step 4 is mainly for the determination of the cut-off value for AACC (MACC). We know that, in the binary treatment case, a well-accepted cut-off value for ASAM is 0.2 [12]. In the continuous treatment case, we set the cut-off value for AACC (MACC) to be 0.1. That is, when AACC<0.1 (MACC<0.1), we claim that the confounding effect between the treatment and the outcome is small. Appendix B shows a heuristic proof for the claimed cut-off value. Figure 1 is an illustration of the Fisher transformation. As can be seen, when |dˉj| is small, the transformation is almost the identity; when |dˉj| is large, |zj| increases faster than |dˉj|. This is another advantage of using Fisher transformation: when we try to minimize AACC, the larger absolute values of correlation coefficients will get more penalty compared to the original scale.

Figure 1

An illustration of Fisher transformation

4.1 Simulation setup

In this section, we conduct simulation studies to compare the performance of the proposed methods to the existing methods. The generation of observed (Y,T,X) is as follows. First, the vector of baseline covariates (potential confounders), denoted as X=(X1,X2,…,X10), is generated from the following distributions: X1,X2∼N(0,1), X3∼Bernoulli(0.5), X4∼Bernoulli(0.3), X5,…,X7∼N(0,1) and X8,…,X10∼Bernoulli(0.5). Among the ten covariates, X1−X4 are real confounders related to both the treatment and the outcome.

The continuous treatment is generated from N(m(X),1), where the mean function is defined for different scenarios. In Scenario (A), m(X) is a linear combination of the real confounders. In Scenario (B), we consider a nonparametric model that is similar to a tree structure with main effects and one quadratic term while in Scenario (C), we add two interaction terms. The true forms of m(X) in different scenarios are displayed as follows:

1. Scenario (A): m(X)=6+0.3X1+0.65X2−0.35X3−0.4X4;

2. Scenario (B): m(X)=6+0.3I{X1>0.5}+0.65I{X2<0}−0.35I{X3=1}−0.4I{X4=0}+0.65I{X1>0}I{X1>1};

3. Scenario (C): m(X)=6+0.3I{X1>0.5}+0.65I{X2<0}−0.35I{X3=1}−0.4I{X4=0}+0.65I{X1>0}I{X1>1}+0.3I{X1>0}I{X4=1}−0.65I{X2>0.3}I{X3=0}.

The potential outcome function for a subject with covariates X is generated from

where ε∼N(0,1). Based on the data generation process, the true dose–response function is

4.2 Results

To compare different methods, we set the value of the parameter of interest to be 0.4, which is the coefficient of T in the dose–response function (11). We apply IPW to estimate the coefficient and employ four different methods to estimate m(X) in the generalized propensity scores: (1) linear approximation (eq. 6) using all the covariates; (2) linear approximation with variable selection; (3) L2 boosting by minimizing the empirical quadratic loss, called mboost by Bühlmann and Yu [22]; (4) boosting with the proposed stopping criteria: Pearson/polyserial, Spearman, Kendall and distance. In Method (2), we employ a variable selection technique that is similar to the idea suggested by Hirano and Imbens [8] to select covariates in the generalized propensity score model. First, we divide the treatment variable into three groups with equal sizes. Then, we test if each covariate is distributed the same in different treatment “groups” using ANOVA at the significance level of 0.05. We only include those covariates that are significantly different among treatment “groups”. In Table 2, we denote Method (1) as Linear¹ and Method (2) as Linear². We generate 1,000 datasets with a sample size of 500. The simulation results are shown in Table 2.

In Scenario (A), the true mean function m(X) is linear in the covariates, and hence the linear approximation proposed by Robins et al. [13] leads to the smallest bias. In addition, Linear² has smaller bias and MSE than Linear¹, which indicates that variable selection in the propensity score model does improve the performance. Compared to Linear¹ and Linear², the proposed methods yield much smaller variances and MSEs, as well as better confidence interval coverages. Compared to the proposed methods, mboost based on L2 boosting yields larger bias and MSE. In Scenarios (B) and (C) where m(X) follows a tree structure, Linear² performs better than Linear¹. In addition, the proposed methods are superior in terms of the bias, MSE and 95% confidence interval coverage. Our simulation results are not very sensitive to the choice of the correlation matrices among Pearson/polyserial, Spearman and Kendall correlations. Distance leads to slightly more biased estimates compared to the other proposed criteria.

To further explore the proposed algorithm, we randomly select 100 datasets from the simulated datasets in each scenario. We then compare the number of trees selected by each criterion with the optimal number of trees that leads to the smallest absolute bias with respect to the true causal effect (0.4). Table 3 shows the average number of trees based on the 100 datasets. Compared to the “best” model with the optimal number of trees, “distance” tends to select smaller models, which explains that “distance” yields relatively larger bias than the other three criteria in the simulation. The differences in the number of trees between the “best” model and the models selected by “Pearson/polyserial”, “Spearman” and “Kendall” are relatively small. In Scenarios (A) and (C), they yield slightly more complex models than the “best” models and in Scenario (B), they yield slightly smaller models.

5.1 Early dieting in girls study

It is reported that dieting increases the likelihood of overeating, weight gain and chronic health problems [23]. We analyze the Early Dieting in Girls study, which is a longitudinal study that aims to examine parental influences on daughters’ growth and development from ages 5 to 15 [24]. The study involves 197 daughters and their mothers, who are from non-Hispanic, White families living in central Pennsylvania. The participants were assessed at five different waves. At each wave, daughters and their mothers were interviewed during a scheduled visit to the laboratory.

In this analysis, we study the influence of mothers’ weight concern on girls’ early dieting behavior. The treatment variable is mother’s overall weight concern (M2WTCON), which is measured at daughter’s age 7. It is the average score of five questions in the questionnaire. A higher value implies the mother is more concerned about gaining weight. In the dataset, its values range from 0 to 3.4. The histogram in Figure 2(A) and the QQ plot in Figure 2(B) show that the treatment is approximately normally distributed. The outcome is whether the daughter diets between ages 7 and 11. We exclude those daughters who reported dieting before age 7, which results in 158 subjects, of which 45 daughters reported dieting. There are 50 potential baseline confounders in this study regarding participants’ characteristics, such as: family history of diabetes and obesity, family income, daughter’s disinhibition, daughter’s body esteem, mother’s perception of mother’s current size and mother’s satisfaction with daughter’s current body [25].

Figure 2

Model diagnostics for the fitted boosting model

5.2 Estimation of the generalized propensity scores

Since the treatment is self-selected, to draw causal inferences, we need to adjust for the confounders in the study. Given a large number of potential confounders, we employ a boosting algorithm to estimate the generalized propensity scores. The simulation studies in Section 4 show that the estimation results are not sensitive to the choice of the correlation matrices among Pearson/polyserial, Kendall and Spearman. Therefore, we use “Pearson/polyserial” as the main criterion to select the optimal number of trees in boosting. Figure 3 displays the AACC value versus the number of trees from 1 to 20,000. Based on the data, the optimal number of trees M=4,846 and AACC=0.11. Figure 2(C) and (D) shows that the residuals from the boosting model (Ti−mˆ(Xi)) have approximately constant variance and they are normally distributed. Based on the residual plots, we conclude that the boosting model sufficiently estimates the treatment-level given covariates.

Now we will focus on assessing the balance in the covariates. Notice that in the original data, AACC=0.177, which is much larger than 0.1, the cut-off value. In addition, if we look at each covariate separately, there are many covariates whose absolute correlation coefficients with T are larger than 0.2. As shown in Figure 4, after applying the weights, most of the absolute correlations among the treatment and each covariate in the weighted sample are below 0.1 on both the original scale and the Fisher transformed scale. This indicates that the confounding effect of the covariates between the treatment and the outcome is greatly reduced after weighting.

Figure 3

AACC value for M=1,…,20,000

Figure 4

Absolute value of correlation coefficients between T and each covariate before weighting (black dots) and after weighting (red plus signs).The left panel shows the original scale and the right panel shows the Fisher transformed scale

5.3 Modeling the mean outcome

To draw causal inferences, the next step is to determine the functional form of the outcome model. The boxplot of the estimated inverse weights from our proposed method is displayed in Figure 5. As shown, most weights are distributed around the value of 1. However, there are some extreme weights with values larger than 10. Extreme weights are harmful to the analysis because they increase the variance of the causal estimates [26]. When we estimate the dose–response function, we shrink the top 5% of the weights to the 95th quantile.

Figure 5

Boxplot of the inverse probability weights using Pearson/polyserial correlation

Since all the potential confounders are well-adjusted by the propensity model, we can model the outcome as a function of the treatment. Otherwise, we may also include covariates that are related to the treatment in the outcome model. We then assume a regression spline function as in eq. (3). For binary outcomes, the regression spline function is

The weighted AIC and BIC displayed in eqs. (4) and (5) are employed to determine the optimal p and K. For a binary outcome, the first part of eqs. (4) and (5) should be replaced by

where pˆi is the estimated probability of early dieting for subject i. We consider three different values for p: p=1,2,3, which corresponds to piecewise linear, quadratic and cubic models. We consider K=0,1,2,…,9. We select the optimal number of p and K based on the values of BICw, which are displayed in Table 4. As shown, the best model is when p=1 and K=0. Therefore, the model we fit is

The causal log odds ratio (β) is estimated as 0.1782 with a standard error 0.2865 (p-value = 0.5349). The standard error is obtained using sandwich formula by the survey package in R. We may also employ bootstrapping to estimate the standard error by repeatedly taking a bootstrap sample with replacement from the original dataset and applying the same estimating procedure. Based on 1,000 replications, the bootstrapping estimate of the standard error is 0.2374, which is slightly smaller than the sandwich formula. It indicates that the probability of daughter’s early dieting increases when the mother’s weight concern increases. However, the causal effect is not significant at the significance level of 0.1.