Discussion of “On Bayesian Estimation of Marginal Structural Models”

Saarela, Stephens, Moodie and Klein (SSMK hereafter) are to be congratulated for writing a really interesting paper. There has been considerable angst on how to reconcile Bayesianity with inverse probability of treatment (IPT) weighting. Fortunately, however, there have been recent insights. For instance, Wang, Parmigiani, and Dominici (2012) provide a model-selection sense in which a treatment-choice model matters for Bayesian inference, while Zigler et al. (2013) shed more light on the thorny issue of “feedback” from the outcome data to the treatment-choice model. And now SSMK may have moved the goalposts considerably. In the subtle longitudinal context, they provide a full and closely argued recipe for Bayesian inference based on modeling the treatment-choice mechanism. To help understand SSMK's method, here it is implemented in a simple example, and then compared to a different method. As a simple technical tool, a Bayesian saturated binary regression model (henceforth BSAT, for short) is defined as follows. For a binary response variable A and binary explanatory variables , a distinct outcome probability is assigned to each of the possible values of B, i.e., . As a conjugate prior, the elements of are taken as iid, with distributions. Hence a posteriori they are independently beta-distributed, and direct Monte Carlo sampling from the posterior distribution is trivial. To bolster understanding of SSMK's proposal, consider the fictitious data in Table 1. These involve subjects, with only two timepoints (), two treatment alternatives (), a binary covariate, and a binary outcome. For ease of exposition, regard () as being “on” (“off”) treatment at the j-th timepoint, while is a harmful outcome, such as having reached irreversible disease progression by the end of the study. As per SSMK, consider the variables as unfolding in the temporal order . In such a small-scale problem there is no need to make parametric modeling assumptions. Each treatment pattern can have a distinct counterfactual mean outcome . Focussing in on the “always-treat” pattern, , the subset of the data pertaining to those treated at both timepoints is summarized in Table 2. Taking these data at face value, without any regard to confounding, yields an estimated outcome probability for the always-take regime to be . Upon fitting BSAT models for , , , and , and then exactly following SSMK's prescription (as per their Supplementary Appendix B), a weight is estimated for each pattern. For these are reported in Table 2, and they yield the weighted data table given in the rightmost columns of the table, i.e., raw cell counts in each row get multiplied by the corresponding weight. The weighted table then gives an estimated outcome probability of for the counterfactual world in which everyone is fully treated. The adjustment for time-varying confounding has quite markedly changed the face-value impression of the data. SSMK's method does not directly use the weighted columns of Table 2, yet reweighting rows of X by Y data tables, for fixed Z, is central. Particularly, the proposal is to Bayesian-bootstrap first, reweight second. Each individual's contribution of exactly one datum is replaced with a contribution of roughly one datum, via the Dirichlet sampling scheme. This gives a “noised up” X by Y data table for fixed Z, which is then reweighted to yield a point estimate of the corresponding . Then the ensemble of such point estimates arising from repeated Bayesian bootstrapping is taken as the posterior distribution of . Doing this, we obtain a posterior mean (SD) of () for . Not surprisingly given the nature of bootstrapping, the posterior mean is effectively the same as the point estimate from the reweighted portion of Table 2. Of course, interest typically lies in comparing different regimes. For , the counterfactual risk difference between fully treating and fully not treating, we get a posterior mean (SD) of (). Stepping back, IPT weighting is but one strategy for such problems. The g-formula approach pioneered by Robins (1986) is another. There are various accounts of this method; see an Appendix to Taubman et al. (2009) for an accessible explanation. The core idea is to probabilistically express the time-evolution of all variables, including both the treatment choice arising under a specific intervention and the treatment choice arising in the absence of any intervention. Assumptions about no unmeasured confounding then equate to assumptions about this joint distribution. In the special case of a single timepoint, the g-formula reduces to the well-known epidemiological procedure of standardization (Snowden, Rose, and Mortimer, ). The key drivers of the scenario are set as , , , and . By dint of all these being positive, the scenario has some typical features. Think of as some manifestation of “being sicker” at the j-th timepoint, as meshes with . Via , those who are sicker are more likely to start treatment, a hallmark of “confounding by indication.” Also, positive and reflect a dual benefit of treatment. There is a direct effect (controlled by ), but also an indirect effect (controlled by ), whereby treatment reduces the chance of the undesirable transition from to . Figure 1 gives results for 50 datasets simulated as described. The previously seen agreement in estimating is no fluke; the SSMK and BGF posterior means are always essentially the same. The corresponding posterior standard deviations agree less closely, though both exhibit very modest variation across repeated sampling. Moreover, in an extended simulation of 1000 datasets, the empirical coverage of 95% equal-tailed credible intervals are 95.6% and 95.7% for SSMK and BGF, respectively, with a discordance rate (one interval covers but the other does not) of only 0.9%. Taking stock then, the two methods start with different premises, and require modeling different parts of the joint distribution of observables. Yet here they give essentially the same estimate and about the same indication of uncertainty. It is comforting to see two seemingly disparate statistical methods, SSMK and BGF, yielding the same answer. In contrasting the methods, however, BGF does have a simple elegance. Models for some conditional distributions must be specified, along with concomitant priors. Then the Bayesian crank is turned, and the target parameter is a function of the unknown parameters. The need for bootstrapping is obviated, and special arguments to support the method's validity are not required.

Fetched live from OpenAlex and de-inverted. Abstracts are not stored in this database: the inverted indexes are 8.6 GB of the frame’s 9.3 GB of text, and the host has 13 GB free.