A Marginal Structural Modeling Approach with Super Learning for a Study on Oral Bisphosphonate Therapy and Atrial Fibrillation

2.2.1 Observed data structure

The observed data on any given patient in this study (referred to as the “BP/AF study” from now on) consist of the collected measurements on exposure, outcome, and confounding variables from study entry until the patient’s end of follow-up. In this analysis, time is measured by consecutive 30-day intervals after study entry. The time when the patient’s longitudinal data are right-censored is denoted by C and is defined by the shortest of (1) the time to the administrative end of study denoted by , (2) the time to disenrollment from the KPNC health plan denoted by , and (3) the time to death denoted by , i.e. . The type of right-censoring event is encoded with a categorical variable denoted by Γ with possible values 1, 2, or 3 to represent end of follow-up by administrative end of study (i.e. ), disenrollment from the health plan (i.e. ), or death (i.e. ), respectively. At each time point , i.e. a 30-day interval poststudy entry, the patient’s exposure to BP is represented by the binary variable , and the patient’s right-censoring status is denoted by the indicator variable . The combination is referred to as the action at time t. At each time point , covariates (e.g. diagnosis of coronary disease or exposure to raloxifene) are denoted by the multidimensional variable and defined as measurements that occur before or are otherwise assumed not to be affected by the actions at time t or thereafter, (). For each time point , the outcome denoted by is defined as the indicator of failure (i.e. a clinically relevant AF event) at time and is an element of the covariates at time . By convention, the outcome at baseline (an element of ) is set to 0, i.e. . To simplify notation, we use overbars to denote covariate and action histories, e.g. a patient’s action history through time t is denoted by . We approached the observed data from the study as realizations of independent and identically distributed (i.i.d.) copies of , where n denotes the sample size. The maximum follow-up time K for any given patient is (i.e. 85 30-day intervals). In Section 2.2.5, we describe how EMR data were concretely mapped into this data structure O for each patient.

2.2.2 Simplified definition of the effect of interest

Informally, a causal effect may be defined by analogy with an ideal randomized experiment, i.e. a trial with perfect compliance and no right-censoring. In this study, the effects of incident one-year exposure to BP on the risk of AF could have ideally been studied by randomizing each patient of the study cohort to one of the following two treatment arms at study entry. In the first arm, patients would have remained unexposed to BP for 1 year, while in the second arm, patients would have been exposed to BP for the same time. The difference between the cumulative risk of AF events at 1 year in each arm is one effect measure that we aim to emulate with observational data in this study.

Formally in the counterfactual statistical framework [24], such causal effects are defined as contrasts between the distributions of potential outcomes in a conceptual experiment. We use notation consistent with that introduced in the previous section to describe the data from this experiment. At each time point of this fictitious experiment, patient’s covariates are observed under any given action intervention denoted by such that right-censoring events do not occur before and at time K, i.e. . These covariates, referred to as counterfactual covariates, are denoted by and satisfy the time ordering assumption formalized by equality: , i.e. actions after cannot affect covariate levels at time t. The counterfactual covariates at , include the counterfactual outcome denoted by . The collection of all counterfactual covariates collected over time during this ideal experiment are referred to as the full data denoted by X, i.e. , where represents the collection of any possible action interventions such that right-censoring events do not occur before and at time K, i.e. . The causal effect that was informally defined earlier by analogy with an ideal randomized trial can now be formally defined as a function of the distribution of these full data:

[1]

where is implied by the full data and represents the counterfactual time to the first failure event under a given action intervention and and represent the two action interventions of interest by which patients do not experience censoring events before and at time 11 (i.e.) while experiencing continuous BP exposure (i.e.) and no BP exposure (i.e. ), respectively. More precisely, if the first failure occurs before time then is defined by the time when the counterfactual outcome jumps to 1 for the first time, i.e. . Otherwise, if failure does not occur during follow-up under action intervention , the counterfactual failure time is defined as . The causal estimand can be viewed as one parameter of a nonparametric MSM for the distribution of the counterfactual survival times .

2.2.3 Estimation approach

Figure 1 represents possible causal relationships between variables collected during the first two periods (on patients followed-up for at least two periods) in the BP/AF study using a causal diagram [25]. Standard estimation methods (e.g. extended Cox regression) to account for confounding involve modeling approaches that compare the distributions of outcomes conditional on both the levels of confounders and exposure variables. It is well-established [26,27,28] that such methods do not properly handle time-dependent confounders affected by early therapy exposure such as covariate L(1) (Figure 1.) because (1) inclusion of such confounders in standard models may lead to bias, since they block some of the causal pathways under investigation; and (2) exclusion of such confounders from standard models may lead to bias due to the lack of confounding adjustment. In observational longitudinal studies that aim to compare the effectiveness of alternate treatment interventions, the existence of time-dependent confounders such as L(1) is often expected since the indications for the treatment decisions under study are also typically risk factors for the clinical outcomes of interest (e.g. CD4 counts and antiretroviral therapy in HIV patients [29] or hemoglobin A1c and treatment intensification in diabetes patients [30]). Thus, standard analytic methods are often suspected to lead to biased inferences in such effectiveness studies. In longitudinal safety studies, however, the existence of time-dependent confounders such as L(1) are often not as obvious since indications for the exposure under study such as history of fracture in the BP/AF study are not typically known risk factors for the clinical outcomes of interest. While time-dependent confounding by covariates such as L(1) is not a priori expected in safety studies, bias from lack of or incorrect adjustment for unforeseen time-dependent confounders remains a concern. For example, it may be alleged that occurrence of major clinical events (e.g. hospitalization for AMI in the BP/AF study) that are risk factors for the outcome under study and also possibly the results of deleterious side effects from prior exposure to treatment will also often lead to the early interruption of that treatment. Such concerns motivate the adoption of an inverse probability of action weighted (IPAW) estimation approach [22, 31, 32] in longitudinal safety studies to not only properly handle possible observed time-dependent confounders on causal pathways of interest but also to adjust for possible time-dependent selection bias due to right-censoring [33].

Figure 1

Causal diagram that represents possible causal relationships between a subset of the variables collected during the first two periods of the BP/AF study (on patients followed-up for at least two periods). T represents the failure time, i.e. time to the first AF event, and represents the indicator of failure occurrence within two periods from study entry. and represent exposure to BP during the first and second periods of follow-up, respectively. represents baseline confounders of the effect of both early, , and late, , BP exposures on failure occurrence within two periods from study entry, . represents intermediate covariates that confound the effect of late BP exposure on failure and that are also affected by the early BP exposure . Standard modeling approaches are not adequate to account for the time-dependent confounding that may result from covariates such as .

Alternatively, one may attempt to alleviate the problem of time-dependent confounding through an intention-to-treat analytic approach where changes in exposure over time are ignored [34, 35]. Confounding adjustment by baseline covariates (L(0)) using standard estimation approaches (e.g. propensity score [36,37,38,39] matching) can then provide unbiased estimates of the effect of the early therapy exposure (A(0)). Such inferences are informative in randomized experiment in which changes in therapies are rare during the follow-up time or else in placebo-controlled trials conducted for evaluating efficacy since inference from an intention to treat (ITT) analysis is then conservative [40]. Placebo-controlled randomized clinical trials conducted for evaluating safety using an ITT approach only can be misleading and have been referred to as “randomized cynical trial” [40], since ITT results in such analyses are anticonservative and may thus conceal deleterious side effects. These considerations support our goal to emulate inference from a placebo-controlled with perfect compliance in the BP/AF study, instead of inference from an ITT analysis of a placebo-controlled trial with possible non-compliance.

Under assumptions detailed below, the causal estimand in the BP/AF study may be estimated consistently based on the following IPAW estimating function^[41]:

[2]

where and g are the generic notation for indicator variables and conditional probabilities of actions, respectively, and where

[3]

represents the earliest time when a right-censoring or failure event occurs or when the action interventions of the hypothetical randomized experiment that we aim to emulate would end (i.e. time point 11). For conciseness, we adopt the following notation

[4]

to represent the difference between the two indicators of whether the patient’s actual exposure to BP and right-censoring status (during the periods that would have been monitored in the ideal trial to be emulated) are concordant with the action interventions in one of the two arms of the ideal trial, i.e. (1) the indicator of both no right-censoring event and continuous exposure to BP up to and at time , and (2) the indicator of both no right-censoring event and continuous non-exposure to BP up to and at time .

As shown in Appendix A, the IPAW estimating function [2] is unbiased at under the following identifiability assumptions [28, 42] which insure that the available observed data contain sufficient information to emulate inference from the ideal randomized trial that would address the research question of interest: (1) the sequential randomization assumption (SRA) formalizes the assumption of no unmeasured confounders in longitudinal studies: . It cannot be tested with data alone. (2) The experimental treatment assignment (ETA) assumption (also known as positivity assumption) insures non-parametric identifiability of the causal effect of interest. In other words, satisfaction of the ETA assumption permits reliable estimation of the causal effect based on true information in the data instead of extrapolation from limited information based on parametric modeling assumptions [43]. Formally, this assumption insures that any patient may experience any of the interventions of interest at any point in time, regardless of her covariate history: and

The following solution of the estimating equations derived from the IPAW estimating function [2] defines an estimator of :

[5]

An estimator of the variance of this IPAW estimator can be derived from its influence curve [24, 44] and can be expressed with the following closed-form formula:

[6]

With observational data, implementation of this IPAW estimator first requires estimation of the product of unknown conditional probabilities of actions g (referred to as action mechanism) and which can be decomposed into two products (referred to as the treatment and right-censoring mechanisms, respectively):

The consistency of the IPAW estimator thus relies not only on the SRA and ETA assumption but also on consistent estimation of the action mechanism, i.e. proper adjustment for observed confounders. The inference drawn from formulas [5] and [6] is asymptotically conservative when the action mechanism is estimated [24]. Bootstrapping can be used instead to draw inference that incorporates variability from estimation of the action mechanism.

Maximum likelihood estimation based on a priori specified and more or less flexible parametric models have been used most commonly in practice to estimate the action mechanism. In few instances [45, 46], data-adaptive estimation of the action mechanism has been implemented as an alternative to the more common reliance on a priori specification of parametric models which often do not reflect true subject-matter knowledge and may thus lead to incorrect inferences [47, 48]. Several machine learning algorithms (also known as learners) are potential candidates for estimating the different components of the action mechanism, e.g. random forests [49], least angle regression [50], logic regression [51], classification and regression trees [52], ridge regression [53], multivariate adaptive regression splines [54], Deletion/Substitution/Addition algorithm [55], high-dimensional propensity score algorithm [56], generalized boosted modeling [57], etc. Akin to the selection of a parametric model, the selection of a machine learning algorithm is rarely based on real subject-matter knowledge about the relative suitability of the different learners available, since “in practice it is generally impossible to know a priori which learner will perform best for a given prediction problem and data set” [58]. A poor choice in machine learning algorithms for estimation of the action mechanism may thus also result in biased inference. Recent theoretical work demonstrated that cross-validation may be used aggressively for estimator selection [59,60,61,62], i.e. to search through a large space of candidate estimators. From this work, a data-adaptive estimation approach referred to as super learning was developed [58]. Its application avoids the need for arbitrary selection of a candidate learner in a particular estimation problem. Instead, super learning combines predicted values from the various candidate learners considered through a weighted average. The selection of the optimal combination of the candidate learners is based on cross-validation to protect against over-fitting such that the resulting learner (called “super learner”) performs asymptotically as well or better than any of the candidate learners considered [58]. Super learning can thus be used for estimation of the action mechanism in this approach to avoid reliance on subjective specification of parametric models.

2.2.4 Extension of the approach over time

The analytic approach just described only makes use of data accumulated during the first 12 follow-up periods. Information in subsequent periods is thus left untapped while it could contribute to addressing the research question of interest. We thus extend the principles of the analytic approach presented earlier to make use of the complete follow-up data available from each patient to evaluate the safety of incident 1-year exposures to BP based on the 1-year cumulative risks of AF.

Informally, the extended approach over time aims to emulate inference from several ideal randomized trials. Each of these trials involves the same two one-year BP interventions described in the previous section but the start of each trial is staggered over time (Figure 2). More specifically, each trial starts at the end of one of the 30-day follow-up period of the AF/BP study and is based on the population of patients who remained unexposed to BP and active in the study (i.e. patients who did not experience a censoring event) through the end of that period. In each of these hypothetical ideal trials, the difference between the cumulative risk of AF events at 1 year in each arm is a measure of the effect of incident 1-year exposure to BP on the risk of AF. The approach described in this section aims to emulate an average of these effect measures by extending the analytic principles described earlier.

To formally describe the generalization of the analytic approach, we extend the overbar and counterfactual notation as follows. The history of a variable (action or covariate) between time t and t′ is denoted by with is nil if . The counterfactual covariate at time t corresponding with an action history that is defined by the following sequence of observed actions and action interventions with such that and is denoted by , where is nil. For patients whose follow-up data are not right-censored through time for (i.e. for whom ), the counterfactual time to first failure at or after t under intervention with any given is denoted by . If failure occurs at or after period t but before period under intervention then is defined by the time when the counterfactual outcome under intervention is 1 for the first time after t, i.e. . Otherwise, if failure does not occur during the interval t to t + s under intervention , the counterfactual failure time is defined as ∞.

The causal estimand and its IPAW estimating function defined by formulas [1] and [2] in the previous section can now be expressed using the new notation above (with and ) as

and

[7]

where . This IPAW estimating functions can be used with observational data from the BP/AF study to emulate inference from a randomized trial that would start at time 0 and end at time 11 with two arms where patients are continuously exposed or continuously not exposed to BP for 12 30-day periods. The inference is generalizable to the target population of patients from which the patients in the BP/AF study cohort were drawn.

The principles of this conceptual ideal trial can now be extended to define another hypothetical ideal trial that would start at time and end at time for any given . The target population for this trial may be defined as the population of patients in the original target population above with the added restriction that patients should not experience BP exposure nor a censoring event before time for any patient i in this new target population). Note that while patients with documented previous exposure to BP are excluded from the target population, patients with a history of AF are included in the target population. The inference from this new trial is also based on the same two action interventions (continuous exposure to BP versus continuous absence of exposure to BP with no censoring event before and at time t denoted by and , respectively). In this new conceptual trial, the causal estimand of interest representing the difference between the cumulative risks of failure in the two arms is defined as:

Inference from such a trial with the available observational data from the BP/AF study may be emulated under the SRA and ETA assumption with IPAW estimation based on the following IPAW estimating function:

[8]

where is defined as the earliest time when a right-censoring or failure event occurs at or after time (start of the trial) or when the action interventions of the trial would end (i.e. time point t):

Note that definition [8] reduces to definition [7] when .

Instead of estimating each effect by solving the IPAW estimating equations associated with separately for each (for a given , we aimed to estimate a single average effect measure ψ by solving the estimating equations associated with

[9]

The solution is the IPAW estimator denoted by of an average measure of the effects of interest for :

[10]

As shown in Appendix B, the average effect measure denoted by that is targeted by this IPAW estimator is given below and can be alternatively defined by a working [63] semi-parametric [64] history-restricted MSM (HRMSM) [65, 66] which, if correct, corresponds with the assumption that the risk differences are constant over time:

[11]

An estimator of the variance of the IPAW estimator _n if the conditional probabilities of actions g (action mechanism) were known can be derived from its influence curve and can be expressed by the following closed-form formula:

[12]

Under the SRA and ETA assumption, the inference for causal estimand from formulas [10] and [12] is asymptotically conservative if the action mechanism is estimated consistently. To avoid reliance on a priori specified parametric models for the action mechanism, super learning may be implemented to select estimators for the different components of the action mechanism as described later. Simulations and applications in previous work [58, 67] have demonstrated the important gains achieved with the super learning approach in prediction problems. To our knowledge, this methodology has however not been implemented in any analyses that approach the scale of observation numbers ( patients followed-up for up to 85 30-day periods) involved in the BP/AF study. Even though we expected a heavy computational burden (computing times and memory requirements) for running the R implementation of the super learner (version 1.1-18) [68] for estimating the action mechanism in the BP/AF study, the reliance on proper confounding adjustment for drawing causal inferences motivated the exploration of the applicability of super learning in a safety study of this scale. We later discuss practical obstacles to the application of super learning in this analysis and practical solutions that permitted preliminary application of its principles.

Figure 2

Visual representation of the principles of the MSM approach in the BP/AF study. The approach aims to emulate inference from ideal randomized trials staggered over time. Each ellipse represents the intervention period for a given hypothetical trial. The green line represents the period during which patients from the original target population for the BP/AF study must remain untreated with BP and active in the study to form the target population of the trial. The red and blue lines represent the action interventions in the two arms of the trial (continuous 1- year exposure to BP is represented in red). The difference between the cumulative AF risks at 1- year in the two arms defines the causal effect measure of interest in each trial and is represented by the arched arrow. The formal definition of this effect using notation from the counterfactual statistical framework is given next to each diagram along with an unbiased IPAW estimating function for that effect.

2.2.5 Implementation challenges and practical solutions

To balance computational burden and misclassification of the time-varying actions and covariates, we chose the 30-day interval as the unit of time in this analysis (see Sections 2.2.1 and 2.2.2). To approximate the observed data structure O resulting from this choice, we first needed to map raw data from the KPNC EMR collected on a continuous time scale to the coarsened time scale chosen for the analysis. The following practical approach was implemented to uphold the time-ordering assumption between measurements of the actions and covariates, i.e. which encodes that covariate represents attribute measurements that occur before the action or that are otherwise assumed not to be affected by the actions at time or thereafter, :

1. Each patient’s follow-up time measured in days from study entry was discretized into consecutive 30-day intervals denoted by . In particular, represents the first days of a patient’s follow-up: and represent the last 30-day interval when a censoring event occurs. Note that unlike t for which represents a monitoring interval of 30 day the follow-up time typically represents a shorter monitoring interval since the actual censoring event often falls within the 30-day interval .

2. The indicator of right-censoring is 0 for and .

3. BP exposure is set to 0 for all time intervals t except for the interval (if any) when BP is first initiated and for all subsequent intervals where BP is used for at least 15 days.

4. For any attribute of the covariate which does not represent a given drug exposure (e.g. diuretic or statins):

   – For all t representing intervals that do not contain the day when BP exposure is first initiated or first discontinued (i.e. when the value for is 1 for the first time or 0 for the first time after a sequence of 1’s, respectively), the value for each covariate attribute of is set to the last measurement of that attribute (if any) (i.e. the last measurement obtained during interval , within the 6 years preceding study entry if the attribute is “cancer.h” and since 1996 for all other attributes.

   – For t representing the interval when BP exposure is first initiated or first discontinued, the value for each covariate attribute of is set to the last measurement of that attribute (if any) or t but always prior to the actual day when BP was initiated/discontinued if , during the first follow-up interval but always prior to the actual day when BP was initiated and otherwise within the 6 years preceding study entry if the attribute is “cancer.h” and since 1996 for all other attributes.

5. For any attribute of the covariates which represents exposure to a given drug (other than BP) at interval t, its measurement was set to 0 except if the drug was used at least 15 days of that interval. This convention relies on the assumption that BP exposure during any given interval should not affect exposure to other drugs during the same interval.

To mitigate computational burden, we implemented a data reduction step based on sampling from the original cohort data. Given the relatively rare onset of treatment with BP during follow-up in the original cohort, i.e.

we differentially sampled a patient’s follow-up data based on whether she was ever treated with BP during follow-up (i.e. ) to insure sufficient information on the treated patients in the resulting reduced data. This complex sampling scheme can be viewed as an artificial censoring event that is superimposed on the original data structure O collected on each patient in the BP/AF study cohort. The resulting reduced data from such complex sampling can be viewed as n i.i.d. copies of , where ∆ is the indicator of a patient’s data being sampled. The probability of a patient’s follow-up data being sampled conditional on her original data is known by design and only based on her “ever/never” exposure status, i.e. , where denotes conditional probabilities of a patient’s data being sampled. We chose the following sampling probabilities by which follow-up data from all patients who ever used BP during follow-up in the original cohort (≈60,000 patients) are retained in the reduced data along with follow-up data from an approximately equal number of patients who never used BP during follow-up in the original cohort: and . The resulting reduced data consist of approximately 8.3 million person-time observations from about 120,000 patients monitored for up to 85 30-day periods each. As shown in Appendix C, the IPAW estimating function for ψ based on the complete data O given by formula [9] can be mapped into the following unbiased estimating function for ψ based on the reduced data O′:

[13]

The resulting reduced-data IPAW estimator of that solves the estimating equations associated with the previous estimating function is

[14]

with . An estimator of the variance of this reduced-data IPAW estimator if the action mechanism g were known can be derived from its influence curve and can be expressed by the following closed-form formula:

[15]

Akin to the estimation function for , estimating functions for the different components of the action mechanism may be mapped into estimating functions based on the reduced data only. In practice, estimation of the action mechanism with the reduced data may thus be implemented using standard estimators (e.g. based on a priori specified parametric models or super learning) where the contribution of each observation is weighted by the inverse probability of sampling to account for sampling bias. Under the SRA and ETA assumption, the inference for causal estimand ψ from formulas [14] and [15] is conservative if the action mechanism is estimated consistently. While bootstrapping can be used instead to draw inference that incorporates variability from estimation of the action mechanism, its implementation was not practical in this study, since it would have required to iterate super learning estimation of each component of the action mechanism on each bootstrap sample.

Estimation of the action mechanism needed as a prerequisite to IPAW estimation in this study was based on further decomposition of both the treatment and right-censoring mechanisms:

Along with the three components of the right-censoring mechanism, only the BP initiation and BP continuation mechanisms need to be estimated for IPAW estimation of the causal estimand . A separate estimator based on pooled data over time j was derived by super learning not only for each of these two components of the treatment mechanism but also for each of the three components of the right-censoring mechanism. Our choice of candidate learners for building these five super learners was mainly based on practical considerations, i.e. the computational burden in terms of both computing speed and memory requirements associated with each learner available in version 1.1-18 of the SuperLearner package in R. The following three types of learners were selected: main term logistic regression,² k-nearest neighbor classification [69], and polychotomous regression [70]. We refer to these estimators as glm, knn, and polyclass learners, respectively, from now on.

For constructing each of the five super learners, six candidate glm learners defined by an increasing number of explanatory variables were considered. More specifically, these glm learners involved the first 13, 20, 30, 40, 50 variables, and all variables, respectively, from an ordered list of explanatory variables. The same list of explanatory variables was considered for building the two super learners predicting BP exposure at time . It consists of the 30-day interval number (named “period”) and all covariates (Tables 1 and 2) associated with interval , excluding the outcome . Note that all explanatory variables are binary variables except for “age” and “period”. The same list of explanatory variables was considered for building the other three super learners and consists of the previous list of covariates but also includes BP exposure at time . For constructing each of the five super learners, the list of relevant explanatory variables was ordered based on not only a combination of subject-matter knowledge but also the magnitude and significance of the crude association between each relevant explanatory variable and the corresponding dependent variable as measured by the odds ratio (OR) and p-value from a univariate logistic regression accounting for repeated measurements.³ More specifically, the ordering of covariates for predicting a given dependent variable was derived based on the following approach: the 14 variables that are established risk factors for AF events (Tables 1 and 2) appear first in the list. They are ordered by decreasing magnitude of association with the dependent variable (as measured by a decreasing distance of the OR from 1). The remaining variables that are significantly associated with the dependent variable at the 0.05 level are ordered by decreasing magnitude of association with the dependent variable and appeared next in the list just before all other variables ordered by decreasing significance of association with the dependent variable (as measured by an increasing p value). Finally, one of the dummy variables for “race” (an established risk factor for AF) is removed from the ordered list due to colinearity with the other dummy variables. Note that the default screening routines available as part of version 1.1-18 of the SuperLearner package in R do not permit the screening of explanatory variables to automatically generate the six glm learners described earlier. Alternate, customized screening routines can however be easily developed and integrated to the super learner algorithm implemented in R.

Due to the heavy computational burden, the polyclass learner considered for constructing each super learner was based on only the 14 variables that are established risk factors for the outcome. We used a customized routine that implements the polyclass learner based on BIC as the model selection criterion to improve computing speed over the default SL.polymar routine that implements the polyclass learner based on cross-validation when the binomial link function is specified.

Given the inability of the default knn learner to incorporate inverse probability of sampling weights, we only considered⁴ the knn learner with k = 10 as a candidate learner when estimating the BP continuation mechanism since the sampling weights for all relevant observation is then 1 by design (i.e. weighting is obsolete). Due to the heavy computational burden, the candidate knn learner was only based on the 14 variables that are established risk factors for the outcome.