Covariate Balancing Inverse Probability Weights for Time-Varying Continuous Interventions

1 Introduction

Propensity Score (PS) methods can be applied to a wide range of different research settings, and have become a popular tool for making causal inference. As shown by Robins [1], PS methods can even be used to evaluate programs that evolve over long periods of time, reacting to changes in the environment that themselves help to bring about. That class of statistical models came to be known as Marginal Structural Models (MSMs) [2], [3], [4]. Despite its theoretical appeal, MSMs as a tool for making causal inference from longitudinal data have yet to extend to general treatment regimes, that is, from binary to (potentially) continuous treatment intensities. Surprisingly, while several extensions of propensity score methods to general treatment regimes have been proposed [5], [2], [6], [7], [8] and keep gaining popularity among applied researchers [9], these methodologies haven’t been coupled with the MSM framework.

Perhaps behind this gap between methodological and applied research is the practical difficulty of estimating the inverse probability weights that are required for MSMs to generate a pseudo-sample from which to estimate the causal quantities of interest. The correct model is generally unknown and its misspecification can actually increase bias, even if the selection on observables assumption holds. In contrast with cross-sectional PS methods, in MSMs the effect of model misspecification is exacerbated, for it propagates across time since the algorithm provided by Robins et al. [2] to estimate the weights in these models involves the product of propensity scores estimated for each time period.

Since, in theory at least, inverse probability weighting renders the intervention independent of time-varying covariates, usual practice involves iteratively checking if weighting on the estimated probability meets this expected property – revising the model until said correlation is minimized. The downside of this algorithm is that it can lead to ad hoc specifications of the probability models which may not agree with the preferred theoretical specification.

As a way to get around this practical difficulty surrounding MSMs, Imai and Ratkovic [10] came up with the idea of formally dovetailing this correlation minimizing property of the weights into the model estimation method itself, eliminating the need to manually and iteratively check the treatment assignment model. They accomplished this by making use of the Generalized Method of Moments (GMM), adding to the model’s score vector – that is, the gradient of the log-likelihood – the correlation minimizing property as an extra moment condition. They have called this new approach the Covariate Balancing Propensity Score (CBPS) [11] and it can easily be implemented through publicly available open-source software [12]. Simulation and empirical studies suggest that the CBPS improves the empirical performance of MSMs by making the treatment assignment model depend less heavily on assumptions which are often not met in practice. However, so far the CBPS extension to MSMs has been confined to binary time-dependent treatments.

Here we propose a further extension to the CBPS methodology that can be used by researchers facing continuous time-dependent treatments. Our proposal builds on the recent generalization of the CBPS methodology to general treatment regimes [13] taking it to the longitudinal setting, thus bridging the gap between two extensions of the original CBPS estimator.

In Section 2 we briefly review MSMs and their assumptions. Section 3 we describe the proposed extension to the CBPS methodology. In Section 4 we present our simulation study. In addition, we present an empirical application in Section 5 to show the performance of the proposed methodology with real data. Finally, we offer concluding comments in Section 6.

2 Marginal structural models

Assume a program that evolves over long periods of time intervening with different intensity across recipients at every stage of their implementation, reacting to changes in the environment which the program itself contributed to bring about; that is to say, reacting to its own previous results in a feedback loop. In this context, variables that can both affect future exposure to the program and be affected by past exposure may confound the impact of the program in a way that no cross-sectional impact evaluation method can sort out [4]. In this Section we briefly review the MSM framework [2], [3], [4] that has been developed to analyze the impact of such programs in case of binary treatments.

Denote recipient i’s observed k-dimensional vector of time-varying covariates at stage t of the program as Xit, with i=1,…,N and t=1,…,T. At every stage t of the program, recipients are observed receiving the exposure level ait to the benefits of the program; that is, one possible realization of the exposure variable Ait, with ai0=0. Collecting all the observed exposures to the program for a given recipient, from its launch up to stage t, gives us the history of exposure a_it=ai1,…,ait. Independent, identically distributed realizations of Yit,A_it,X_it are observed for each recipient at every stage of the program, where Yit is some outcome variable of interest. Let X_it and x_it be similarly defined for a covariate history, where Xit is the most recent set of variables that could affect Ait, and are also affected by past exposure history, A_it−1.

Following Robins’ adaptation of the potential outcome framework [14], [15] to longitudinal settings [16], we use Yita_t to represent the potential outcome for recipient i measured at time t associated to the exposure history a_ up to t, which we assume does not depend on the exposure histories of other recipients. This is often referred to in the literature as the stable unit treatment value assumption. Naturally, any recipient of the program exhibits only one of these potential outcomes at any given stage t; that is, the one associated to its own particular exposure history up to that point in time Yit=Yita_it, the consistency assumption. All other, unobserved, potential outcomes are said to be counterfactual.

As with other PS methods, MSMs are based on the assumption of no unmeasured confounding, in this case at each time period. This assumption, also referred to as sequential ignorability [4] or conditional interchangeability [3], states that the exposure level at stage t is statistically independent (ignorable) of the potential outcomes, conditional on the covariate and exposure histories up to that point. Using previous notation, sequential ignorability says that for any exposure history up to t, at, Yita_t⊥Ait∣A_it−1=a_t−1,X_it=x_t. Additionally, also typical in PS methods, MSMs also assume that at any stage t, there is no covariate history x_t and past exposure a_t−1 such that all recipients with such histories are certain to receive the identical exposure at. That is, each exposure history must have some positive probability of occurring. This last assumption can be formally written as, if fA_it−1,X_it≠0 then fA_it∣A_it−1,X_it>0, where f stands for the probability density function (pdf).

Under the above assumptions (sequential ignorability, consistency and positivity), Robins [1] showed that the conditional probability density of exposure history can be used to consistently estimate the impact of such dynamic interventions. In this generalization of Inverse Probability of Treatment Weights (IPTW), the probability of observing a particular exposure history is expressed as the product of the respective probability density at every stage, ∏t−1TfAt∣A_t−1,X_t. However, using this probability directly as weights when fitting the outcome model leads to highly variable estimates and researchers generally follow Cole and Hernán [17] and use stabilized weights of the form

where, for convenience, we define SW0=1.

In nonrandomized evaluations, irrespective of the nature of the intervention whose impact one might be trying to estimate – whether it’s best coded as binary, multivalued or continuous –, these weights are unknown, and have to be estimated. Usually, a parametric model is used to estimate the probability density of exposure given the past at each time period. When the intervention comes as a single action that either occurs or not at every stage of the program (i. e. coded as a binary variable), a common approach is to estimate this probability density at every stage with a standard logistic regression (logit) model. In the case of continuous interventions, Gaussian models (with constant variance) are often assumed, i. e., At|A_t−1,X_t∼NhA_t−1,X_t;γ,σ2, which can be estimated by the GMM.^[1] If there is reason to believe that h is a linear, additive function of exposure and covariates histories, and parameters θt=α,β,σ2, researchers may regress the current level of exposure on all past exposure and covariates,^[2] e. g.,

The numerator of the stabilized weights is typically estimated using the empirical distribution of the exposure. If the assumptions behind equation 2 hold, and further assuming the marginal distribution of At to be Gaussian with mean zero (possibly after the pertinent transformation), the stabilizing weight is given by^[3]

Once these stabilized weights are estimated, all the researchers have to do to is use them in running whatever model of the outcome, as function of treatment history alone, they would have used in the case of a randomized evaluation.^[4]

3 Correlation breaking IPTW

As with other PS methods, Robins’ stabilized weights break the link between the exposure history and the time-varying covariates, guaranteeing the exposure history is unconfounded. That is, at any given stage of the program, we expect the histories of both exposure and the time-varying covariates to be uncorrelated in the weighted data [4], [18], [19].

Formally, under our time-varying setting, we can see this correlation-breaking property for all t as

Note that if we center both the exposure variable and the covariates, equation 4 equals zero. We propose to use this weighted cross moment as an extra moment condition in a GMM framework, which for t=1 is the same Fong et al. [13] used in the cross-sectional setting.

Including equation 4 among the score conditions – the gradient of the log-likelihood with respect to θ – as moments, allows to account for both properties we would expect the stabilized weights satisfy in a single estimator: predict the exposure among recipients and break the link between the exposure and the time-varying covariates.

4 Simulation studies

In order to assess the empirical performance of these correlation-breaking stabilized weights, we carry out some simulation studies. In particular, we consider the case of three time periods, i. e., t=1,2,3, and use five different sample sizes n= 300, 600, 900, 1200 and 1500, corresponding to 100, 200, 300, 400 and 500 exposure-recipients across our time frame. In our setup, summarized in Figure 1, the treatment-generating process is a function of a vector of exogenous time-varying covariates, X_it, and the previous history of exposure, A_it−1. The outcome variable, Y, is a function of both entire histories of covariates and exposure.

Figure 1

Data Generating Process in Simulation Study.

We use this data generating process (DGP) to examine the performance of the proposed methodology, vis-à-vis the usual Ordinary Least Squares (OLS) approach, when faced with both misspecification of the treatment-assignment model – incorrectly specifying the lag structure, distributional family and cross-time homoscedasticity – and measurement error.

We use a similar simulation setup as Imai and Ratkovic [10] in two scenarios. In the first scenario, the simulation consists of four covariates. For time t, we use the covariates

where each Zitk is an i.i.d. draw from the standard normal distribution, and Uit is constructed as

for t=2,3 and Uit=1 for t=1. The continuous treatment Ait is generated by

where α=.8,−.2,.4,−.6⊤ and vi is a standard normal disturbance – a cross-time homoscedastic normal model. All parameters were selected in an almost random manner looking only to guarantee a minimum dispersion. The continuous outcome Y is generated by a linear combination of the whole exposure and confounders histories such that

where δ=−23.5,12.5,7.3,−17.4⊤ and ui is another normal disturbance with mean zero and standard deviation five. Note that we have generated the outcome variable for every t, making it a function of both entire histories of covariates and exposure.^[8]

There are many ways in which researchers can misspecify the treatment assignment model. In our simulations we explore two scenarios. In the first misspecified scenario, we use the following non-linear transformation of the covariates

as measurement error, as well as incorrectly specify the lag structure of the treatment-assignment model using only the data from the immediately previous period such that

but maintain the correct outcome model. We report results of 1,000 simulated datasets.

In this first scenario, with every dataset, we fit a Pooled OLS model (GLM) as the treatment assignment model, using correct and incorrect model specifications as discussed above. Then we fit the same models using the proposed correlation-breaking IPTW methodology (CBIPTW) in two ways: first using both score and correlation-breaking conditions (Over-identified Correlation Breaking Inverse Probability of Treatment Weights (CBIPTW)) as described in section 2, and with the correlation-breaking condition alone (Just-identified CBIPTW), ignoring the score conditions. This gives us three different estimates according to the importance placed on both expected properties of the weights. The Generalized Linear Model (GLM) estimates obviate the correlation-breaking property, ignoring the subvector g3, while the just-identified CBIPTW pays no attention to predict exposure assignment, ignoring the subvector g1,g2⊤. Over-identified CBIPTW stands in the middle attending to both criteria.^[9] Finally, we regress in a pooled model the outcome variable on the cumulated exposure, ∑t′=1tAit′ using the stabilized MSM weights. The coefficient that results from this last regression is compared with the numerical estimate obtained from a thousand simulations using large datasets (n=30,000) and the true exposure assignment probabilities.

Table 1 presents the results for our five sample sizes. The columns show the bias and Root Mean Squared Error (RMSE) from the pooled weighted linear regression of the outcome on the cumulated exposure using the stabilized weights. The results are in line with those of previous Covariate Balance Propensity Scores for MSMs [10]. When the exposure assignment model is correctly specified, all methods have small bias and small RMSE. In contrast, when the model is misspecified, the bias and RMSEs are large and even grow with the sample size. However, the RMSE of the CBIPTW estimates grow at a much smaller rate, thereby outperforming OLS.

In the second misspecified scenario, we use almost the same DGP, only this time vi is a Gamma disturbance with different variance at every time-point: .5, 1 and 2 for t1, t2 and t3 respectively. In this case, the misspecification consists in fitting a normal cross-time homoscedastic model for the treatment Ait, with no measurement error and the correct lag structure.

Table 2 shows the result of this second misspecified scenario, fitting a Gaussian cross-time homoscedastic treatment-assignment model, just as in the first scenario, but this time to a Gamma, cross-time heteroscedastic DGP. Here too we see the RMSEs of the CBIPTW grow at a much smaller rate, outperforming again the traditional OLS approach, this time also in bias.

What the statistical evidence shows is that the proposed method effectively improves the data fitting processes narrowing the estimator’s dispersion, albeit at the expense of bias, resulting in more precise estimates than GLM. In a way, the extra balancing conditions work as a muffling effect that prevents the OLS estimates (that only consider the score conditions) from drifting away whenever the model is misspecified.

In short, Tables 1 and 2 tell a simple story for time-varying continuous interventions: if researchers are lucky enough to know the true DGP, and it happens to be Normal, they are probably better off with OLS. If, on the other hand, researchers have reasonable doubts with respect to the origin of their data, they are probably better off with CBIPTW, given that it is less affected by (small) departures from the assumptions underlying the traditional OLS approach.

Together, what these two sets of simulations show is that the balancing conditions counterweight quite general misspecifications of the treatment-assignment probability model at roughly the same bias – or less – and quickly stabilizing an otherwise exploding RMSE.

5 Empirical application

As an empirical application of the proposed method, we estimate the impact of the Mexican universal health insurance program, Seguro Popular (SP), on key variables associated with the provision of healthcare services by the Secretaría de Salud – Ministry of Health – (SSA). We analyze data from the federal records on infrastructure and human resources employed by the SSA,^[10] and the coverage records of the program at a yearly frequency from 2001, the year prior the introduction of the program, to 2013. The information is aggregated to the Sanitary Jurisdiction (SJ) level, the most basic administrative units of SSA in charge of the operation of healthcare services and its programs. Time-varying covariates include doctors’ offices, staffed and non-staffed hospital beds, physicians and nurses, both with and without day-to-day contact with patients, all per 1,000 people outside the social security network.^[11] In order to make the common support assumption more credible and the estimates less model dependent, 39 out of 233 SJs are dropped from the analysis.^[12]

Figure 2

Covariate-Balance for Time-varying SP coverage on Sanitary Jurisdictions.

Figure 2

(continued).

The estimations of the treatment exposure model fit a single probability model, regressing the yearly increment of SP coverage (per 1,000 people outside the social security network) on the most recent set of covariates and past (cumulative) SP coverage.^[13] Looking at the degree of covariate balance achieved by the CBIPTW (Fig. 2, Appendix Table A.1), we see that it is improved in all cases, as medians are closer to zero, and that is best achieved using the just-identified approach, both on location and dispersion. It is worth noting that, in this case, fitting only the probability model (GLM) does not break the correlation between the exposure to treatment and three covariates: doctors’ offices, physicians and nurses in day-to-day contact with patients.

Using weighted linear regressions to estimate the impact of (cumulated) SP coverage on the change in the (relative) numbers of doctors’ offices, physicians and nurses (per 1,000 people outside the social security network) providing clinical care, we can see (Fig. 3, Appendix Table A.2) that in all probability the effect of SP is positive in all three outcome variables, and grows bigger the less correlated are the covariates with the treatment in the pseudo-sample.

Analysis by Bosch and Campos-Vzquez [28] on the impact of SP on the expansion of healthcare services ignored both the duration of the roll-out process and the different levels of treatment intensities observed during that process, applying a differences-in-differences approach where treatment was dichotomized (municipalities were considered as ‘treated’ once more than 10 individuals had been affiliated to SP). Huffman and van Gameren [19] improve on that by both considering a continuous treatment indicator (as in this paper) as well as a better-suited unit of analysis. While our results are somewhat similar to earlier findings, to the extent that positive impacts of SP on healthcare services is reported, our results – slightly bigger than Huffman and van Gameren [19] – should be considered more precise as we use an estimator both less model dependent and that outperforms the balance achieved trough GLM.

It is important to note that, in empirical applications, the least comparable the units are, that is to say, the least reasonable the common support assumption is, the more relevant the CBIPTW become. This is because the CBIPTW automatically compensate for the instability in the probability model generated by the lack of common support, something King and Zeng [24] identified as the risk of extreme counterfactuals. In this sense the CBIPTW may also contribute to expand the validity of the results, allowing the analysis of larger data sets that otherwise would be unreasonable.

Figure 3

Average Treatment Effects of SP coverage on Sanitary Jurisdictions.

6 Conclusions

We have shown that the idea behind the Covariate Balancing Propensity Score [11] is flexible enough to cover also longitudinal settings with general treatment regimes. Here we provide an estimation method for Marginal Structural Models with continuous treatments whose main advantage is that it directly includes the correlation-breaking property of Robins’ [1] stabilized weights in the estimation procedure. This new tool avoids the manual process of checking the correlation amongst the current exposure to treatment and the covariate histories in the weighted data and then respecifying. This new extension to the CBPS may help clear the way for those interested in causal analysis in these complex settings.^[14]

Since the balance conditions do not require making any modeling assumptions in themselves, but only that the weights balance the covariates – something we know to be true regardless of the true DGP –, it can be argued that minimizing the correlation between covariates can be done more efficiently then with already available methods, such as the non-parametric extension of the CBPS to general treatment regimes [12] or Genetic Matching [29]. This point grows in importance as treatment assignment models gain in complexity.^[15] However, unlike these other methods, the methodology we propose here provides the researcher also with a probability model that outperforms the traditional GLM approach from which it is possible to start building bridges between MSMs and the decomposition methods literature, something that remains a pending task for future research.