Estimating the causal effects of multiple intermittent treatments with application to COVID-19.

To draw real-world evidence about the comparative effectiveness of multiple time-varying treatment regimens on patient survival, we develop a joint marginal structural proportional hazards model and novel weighting schemes in continuous time to account for time-varying confounding and censoring. Our methods formulate complex longitudinal treatments with multiple ``start/stop'' switches as the recurrent events with discontinuous intervals of treatment eligibility. We derive the weights in continuous time to handle a complex longitudinal dataset on its own terms, without the need to discretize or artificially align the measurement times. We further propose using machine learning models designed for censored survival data with time-varying covariates and the kernel function estimator of the baseline intensity to efficiently estimate the continuous-time weights. Our simulations demonstrate that the proposed methods provide better bias reduction and nominal coverage probability when analyzing observational longitudinal survival data with irregularly spaced time intervals, compared to conventional methods that require aligned measurement time points. We apply the proposed methods to a large-scale COVID-19 dataset to estimate the causal effects of several COVID-19 treatment strategies on in-hospital mortality or ICU admission, and provide new insights relative to findings from randomized trials.

Introduction

The COVID-19 pandemic has been a rapidly evolving crisis challenging global health and economies. Public health experts believe that this pandemic has no true precedent in modern times. While multiple COVID-19 vaccines have been developed across the globe, no consensus has been reached on optimal clinical management of COVID-19 (Yousefi et al., 2020). The lack of evidence for effective treatment options warrants further investigation into the causal effects of multiple COVID-19 treatment strategies currently implemented in clinics. Although randomized controlled trials (RCTs) are considered as the gold standard for evaluating the efficacy of COVID-19 therapies, they are enormously expensive and time consuming, especially in a time of crisis. Stringent inclusion and exclusion criteria also limit the generalizability of RCTs to frailer populations at higher risk for severe morbidity and mortality. To overcome these challenges, we study the causal effects of COVID-19 treatment strategies on patient survival by leveraging the continuously growing observational data collected at the Mount Sinai Health System—New York City’s largest academic medical system. We focus on four commonly used medication classes that are of most clinical interest: (i) remdesivir; (ii) dexamethasone; (iii) anti-inflammatory medications other than corticosteroids; and (iv) corticosteroids other than dexamethasone.

The complex nature of COVID-19 treatments, owing to differential physician preferences and variability of treatment choices attributable to evolving clinical guidelines, poses three major challenges for statistical analysis of observational data that cannot be easily addressed by existing longitudinal causal inference methods. First, treatment is not randomly allocated and the treatment status over time may depend upon the evolving patient- and disease-specific covariates. Second, the measurement time points during the follow-up are irregularly spaced. Third, there is more than one treatment under consideration. Patients can be simultaneously prescribed to various treatment combinations, or can be switched to a different treatment. Figure 1 illustrates the observed treatment trajectories for nine randomly selected patients during their hospital stays.

Fig. 1.

Treatment processes for nine randomly selected patients visualized by heat maps. Colors indicate remaining on treatment. Lack of color corresponds to being switched off treatment.

While previous work has shown that a continuous-time marginal structural model is effective in addressing time-varying confounding and provides consistent causal effect estimators (Johnson and Tsiatis, 2005; Saarela and Liu, 2016; Hu et al., 2018; Hu and Hogan, 2019; Ryalen et al., 2020), the development has been restricted to a single longitudinal treatment and therefore may not be directly applicable. We consider a joint marginal structural model to accommodate multiple longitudinal treatments in continuous time. To estimate causal parameters in the joint marginal structural model, we derive a novel set of continuous-time stabilized inverse probability weights by casting each treatment process as a counting process for recurrent events, allowing for discontinuous intervals of eligibility. In addition, we propose to use machine learning and smoothing techniques designed for censored survival data to estimate such complex weights. Through simulations, we demonstrate that our approach provides valid causal effect estimates and can considerably alleviate the consequence of unstable inverse probability weights under parametric formulations. We further undertake a detailed analysis of a large longitudinal registry data of clinical management and patient outcomes to investigate the comparative effectiveness of multiple COVID-19 treatments on patient survival.

Joint Marginal Structural Survival Model

Notation and set up

We consider a longitudinal observational study with multiple treatments and a right-censored survival outcome. Denote t as the time elapsed from study entry (e.g., hospital admission), t the maximum follow-up time, and a collection of time points on the interval [0, t]. Suppose each individual has a p-dimensional covariate process , some elements of which may be time-varying; by definition, the time-fixed elements of L(t) are constant over . Let T denote time to an outcome event of interest such as death, with as its associated zero-one counting process. We consider W different medication classes (treatments), whose separate and joint causal effects on patient survival are of interest. We use A(t) to denote the assignment of treatment , which can be characterised as the counting process, where we let A(t) = 1 if an individual is treated with w at time t and A(t) = 0 otherwise (Johnson and Tsiatis, 2005). Let C denote the time to censoring due to, for example, discharge or loss to follow up. We use the overbar notation to represent the history of a random variable, for example, corresponds to the history of treatment A from hospital admission up to time t and corresponds to the covariate history up to time t. Following the convention in the longitudinal causal inference literature (Robins et al., 2008), we assume the treatment decision is made only after observing the most recent covariate information just prior to the treatment; that is, for a given t, A(t) occurs after L(t) for all w.

Let represent the counterfactual failure time to event of interest had an individual been assigned treatment history rather than the observed treatment history . Similarly, represents the observed failure time to event for an individual given the observed treatment history. We similarly define as the counterfactual censoring time under treatment . The observed data available for drawing inferences about the distribution of potential outcomes are as follows: the observed time to outcome event is T* = T ∧ C, with the censoring indicator Δ = I(T ⩽ C). Note that both the treatment processes {A(t), w = 1, …, W} and the covariate process are defined for all but are observed only at discrete and potentially irregularly spaced time points for each individual. For example, individual i may have covariates and treatment status observed at a set of discrete time points from study entry t = 0 to his or her last follow-up time . We denote the set of discrete time points with observed covariate and treatment information for individual i as , and therefore the observed covariate and treatment histories become and .

Joint marginal structural model for survival outcomes

We consider a marginal structural model to estimate the joint causal effects of on patient survival. The most popular model specification is a marginal structural Cox model, for its flexibility in handling baseline hazard and straightforward software implementation when used in conjunction with the stabilized inverse probability weights (Howe et al., 2012). When there is a strong concern that the proportional hazards assumption may not be satisfied across the marginal distribution of the counterfactual survival times, alternative strategies including the structural additive hazards model or accelerated failure time model can also be considered. For purposes of presenting our methodology, we focus on the marginal structural Cox model but extensions to alternative structural models are possible. For notational brevity but without loss of generality, we first consider W = 2 treatments. Expansion of the joint marginal structural model and weighting schema for W ⩾ 3 treatments is discussed in Section 3.4. Specifically, we assume follows a marginal structural proportional hazards model of the form where is the hazard function for and λ0(t) is the unspecified baseline hazard function when treatment A1 and A2 are withheld during the study. The parameter ψ1 encodes the instantaneous effect of treatment A1 on in terms of log hazard ratio while A2 is withheld during the study. Similarly, ψ2 corresponds to the instantaneous treatment effect for A2 in the absence of A1. The multiplicative interaction effect of A1 and A2 is captured by ψ3. In addition, the hazard function can depend on baseline covariates by elaborating model (2.1) or by using a stratified version of λ0(t). Model (2.1) implicitly assumes that the instantaneous treatment effect is constant in the course of follow-up. This model assumption is reasonable given that the COVID-related hospitalization is generally short and medications are prescribed for days in succession. Finally, model (2.1) is a continuous-time generalization of the discrete-time model considered by Howe et al. (2012) for estimating the joint survival effects of multiple time-varying treatments.

Model (2.1) offers two advantages for the estimation of treatment effects. First, the counterfactual survival function can be expressed as . Therefore, causal contrasts can be performed based on any relevant summary measures of the counterfactual survival curves including median survival times and restricted mean survival times. Second, model (2.1) allows for the estimation of causal effects of interventions defined by varying treatment initiation timing and treatment duration. For example, an intervention may take the form of , representing prescribing treatment A1 until t† (e.g., t† = day 6). A more complex intervention strategy is , which refers to assigning treatment A1 until t†/2 and then switching altogether to A2 until t†.

Estimating Structural Model Parameters in Continuous Time

To obtain a consistent estimator for ψ = {ψ1, ψ2, ψ3} in model (2.1) using longitudinal observational data with two treatments, we introduce the following causal assumptions and maintain them throughout the rest of the article:

(A1) Consistency. The observed failure times, , where . Similarly for the observed censoring times, . The consistency assumption implies that the observed outcome corresponds to the counterfactual outcome under a specific joint treatment trajectory when an individual actually follows treatment . This is an extension of the consistency assumption developed with a single time-varying treatment (Robins, 1999) to two time-varying treatments.

(A2) Conditional Exchangeability. Alternatively referred to as sequential randomization, this assumption states that initiation of treatment at time t among those who are still alive and remain in the study is conditionally independent of the counterfactual survival time conditional on observed treatment and covariate histories. Mathematically, let denote the observed history up to t−, then where is the joint intensity process of the joint counting process A1(t) and A2(t). Similarly, let denote the observed history up to t, we assume conditional exchangeability for censoring such that , where λ(t) is the intensity process corresponding to the counting process of censoring. Our conditional exchangeability assumption is a continuous-time generalization of the usual sequential randomization assumption for the discrete-time marginal structural models (Robins, 1999; Howe et al., 2012).

(A3) Positivity. We assume that at any given time t, there is a positive probability of initiating a treatment plan, among those who are subject to initiating at least one treatment, for all configurations . For a pair of joint treatments (A1, A2), at a given time t, individuals with treatment status (0, 0), (0, 1) or (1, 0) are subject to “initating” at least one treatment. The treatment initiation patterns can be as follows: (0, 0) → {(0, 1), (1, 0), (1, 1)}, (0, 1) → {(1, 1),(1, 0)}, (1, 0) → {(1, 1), (0, 1)}. Treatment discontinuation (from 1 to 0), however, is not considered a stochastic process in our study, as COVID medication is typically prescribed with a specific treatment duration, e.g., treat with dexamethasone at a dose of 6 mg once daily for 10 days (RECOVERY Collaborative Group, 2021). Furthermore, because an individual cannot be at risk for receiving the same treatment once he or she is on the treatment, we only need to assume the positivity when the individual is off that specific treatment, i.e., at risk for initiating that treatment.

Framing repeated treatment initiation as recurrent events

As Figure 1 suggests, the observed treatment pattern is complex due to considerable variability in COVID treatment protocols and clinician preferences over time. Individuals may discontinue a treatment and restart the same treatment at a later time; or they may be switched altogether to another treatment. Meanwhile, patients can take more than one treatment for a period of time. Each treatment can therefore be viewed as the counting process of recurrent events, with discontinuous intervals of treatment eligibility (Andersen and Gill, 1982). Specifically, casting treatment initiation as a recurrent event process captures two distinguishing features of our observational data: (i) having received a treatment would prevent an individual from receiving the same treatment again for the time period while the individual is on the treatment; and (ii) after the individual was off the treatment, he or she would be eligible or at risk for re-initiating the treatment.

To formalize the treatment initiation process, we first consider a univariate treatment process . We assume that the jumps of A(t), i.e., dA(t), is observed on certain subintervals of only. Specifically for individual i, we observe the stochastic process A(t) on a set of intervals , where . This representation implies the following results. First, an individual can have at most J ⩾ 1 initiations of treatment w: if , then individual i has J − 1 treatment initiations; and if , then individual i has J treatment initiations. A special case where J = 1 and corresponds to the situation where individual i is continuously eligible for treatment initiation and has not been treated with w during the follow-up. Second, once treatment is initiated, A(t) is no longer stochastic until person i discontinues the treatment. This also suggests that the jth treatment initiation is observed at U. Third, A(t) = 1 ∀ t ∈ (U, V], j = 1, …, J − 1. In words, treatment status is equal to one deterministically on the discontinuous intervals of ineligibility (i.e., on treatment period). Define a censoring or filtering process by , and the filtered counting process by , where q indexes the qth treatment initiation.

Following Andersen et al. (1993), we assume conditional independence among occurrences of treatment initiation given all observed history, and that the set is defined such that is predictable. The observed data with occurrences on the set can therefore be viewed as a marked point process generating the filtration . Similarly, we denote the filtration generated by the counting process corresponding to by . We assume A(t) follows Aalen’s multiplicative intensity model (Aalen et al., 2008) λ(t, θ) = α(θ, t)Y(t), with respect to , where λ(t, θ) is the intensity process of A(t), α(θ, t) is the hazard rate function parameterized by θ, and Y(t) is the at-risk function with Y(t) = 1 indicating person i is eligible just before time t for the qth initiation of treatment w in the interval [t, t + dt), and Y(t) = 0 indicating otherwise. It follows that the filtered counting process follows the multiplicative intensity model with respect to (Andersen et al., 1993). Here, . With two treatments, model (3.4) can be directly extended for the joint treatment initiation process as where is the at-risk process for the qth treatment initiation with the filtering process defined jointly by A1 and A2.

Derivation of the continuous-time weights

Under assumptions (A1)–(A3), a consistent estimator of ψ can be obtained by solving the weighted partial score equations (Hu et al., 2018), where is the weight that corrects for potential time-varying confounding for time-varying treatments A1 and A2, Z(A1, A2, t)(3×1) = [A1(t), A2(t), A1(t)A2(t)]⊤, and is a modified version of the weighted mean of Z over observations still at risk for the outcome event at time t. In equation (3.7), we define the weighted risk set indicator for outcome , where is the at-risk function for the outcome event, and r(a1, a2, t) = exp{ψ1a1(t) + ψ2(t)a2(t) + ψ3a1(t)a2(t)}.

In the discrete-time setting with non-recurrent treatment initiation, the stabilized inverse probability weights (we suppress subscript i for brevity) are given by Howe et al. (2012) where t’s are a set of ordered discrete time points common to all individuals satisfying 0 = t0 < t1 < t2 < … ⩽ t. While in (3.8) corrects for time-varying confounding by adjusting for in the weights, it requires that the time points are well aligned across all individuals. In addition, it does not accommodate the recurrent nature of complex intervention strategies as in our observational study.

We now generalize the weights developed for the discrete-time setting to a continuous-time process, which do not require the time points to be well aligned. Partition the time interval [0, t] into a number of small time intervals, and let dA(s) be the increment of A over the small time interval [s, s + ds), ∀s ∈ [0, t]. Recall that treatment initiation, or the jumps of A(t), dA(t), is observed on a number of subintervals of only. That is, conditional on history , the occurrence of treatment initiation for an individual in is a Bernoulli trial with outcomes dA(s) = 1 and dA(s) = 0. Then the P (A(t) = a(t)|•) in equation (3.8) can be represented by which takes the form of the individual partial likelihood for the filtered counting process . When the number of time intervals in [0, t] increases and ds approaches zero, the finite product over the number of time intervals of the individual partial likelihood will approach a product integral (Aalen et al., 2008), given by where ΔA(t) = A(t)−A(t−). For individual i, both factors in (3.9) need to be evaluated with respect to the individual’s filtered counting process , with the first quantity being equal to the finite product over the jump times and the second quantity being the survival function for treatment initiation. As described in Section 3.1, the number of treatment initiations for individual i, Q can take three values: (i) Q = 0, (ii) Q = J − 1 or (iii) Q = J. Corresponding to the three cases, the quantity in (3.9) can be rewritten as where and are the survival and density function of the filtered counting process for treatment A. Here we assume that initiations of different treatments are ordered. For example, whether to initiate A2 at time t is decided upon observing the treatment status A1(t). This suggests that the hazard function is estimable by conditioning on and ; and the hazard function is estimable by conditioning on and . For exposition brevity, let , , and . Putting this all together, the individual continuous-time stabilized inverse probability weight that corrects for time-varying confounding by is given by with being:

Turning to censoring, under the conditional exchangeability assumption (A2), the censoring process is covariate- and treatment-dependent. To correct for selection bias due to censoring, we additionally define a weight function associated with censoring, where S is the survival function associated with the censoring process, and This leads to a final modification of the estimating equation for ψ, where and .

Estimation of the causal survival effects

We consider four ways in which the continuous-time weights can be estimated: (i) fitting a usual Cox regression model for the intensity process of the counting process of treatment initiation , estimating the density function and survival function from the fitted model with the Nelson-Aalen estimator for the baseline intensity function, and finally calculating the weights following equation (3.10); (ii) smoothing the Nelson-Aalen estimator and in turn and estimated from the fitted Cox regression model by means of kernel functions (Ramlau-Hansen, 1983), and calculating the weights using the smoothed version of and ; (iii) fitting a multiplicative intensity tree-based model (Yao et al., 2020) in which the functional form of the intensity ratio for treatment initiation is flexibly captured to estimate the and and calculate the weights; (iv) smoothing the Nelson-Aalen estimator of the baseline intensity from the tree-based model and calculating the weights using the smoothed version of and . Among these approaches, (i) relies on the parametric assumptions about the intensity ratio relationships between the treatment initiation process and covariate process and may be subject to model misspecification and bias for estimating causal effects. Compared to the Nelson-Aalen estimator which includes discrete jumps at event occurrences, the kernel function estimator in (ii) may help alleviate the issue of extreme or spiky weights, and has also been shown to be a consistent and asymptotically normal baseline intensity estimator (Andersen et al., 1993). Approach (iii) leverages a recent random survival forests model (Yao et al., 2020) that can accommodate time-varying covariates to mitigate the parametric assumptions and attendant biases associated with the usual Cox regression. With baseline time-fixed treatment, prior work has used similar machine learning techniques to improve propensity score weighting estimators (Lee et al., 2010) as well as to provide more accurate causal effect estimates with censored survival data (Hu et al., 2021). Finally, approach (iv) smooths the baseline intensity estimated from the survival forests for estimating the stabilized inverse probability weights, and serves as an additional step to smooth over the potentially spiky weights. In Section 4, we compare the performances of these four strategies to estimating the continuous-time weights to generate practical recommendations. In addition, the censoring weight function Ω(G) can be estimated in a similar fashion via any one of these four approaches. Additional details of kernel function smoothing in approach (ii) and random survival forests in approach (iii) are presented in Web Appendix S1.

To accommodate the time-varying covariate process and account for the recurrent nature of treatment initiation, we fit a survival model to the counting process style of data input. Each individual is represented by several rows of data corresponding to nonoverlapping time intervals of the form (start, stop]. To allow for discontinuous intervals of eligibility, when defining multiple time intervals on for individual i, the duration of a treatment is removed from when the individual is currently being treated and therefore no longer eligible for initiating the treatment. Finally, since our estimators for ψ is a solution to the weighted partial score equation (3.11), we can use the robust sandwich variance estimator to construct confidence intervals for the structural parameters. In practice, the robust sandwich variance estimator is at most conservative under the discrete-time setting (Shu et al., 2021), and we will empirically assess the accuracy of this variance estimator with continuous-time weights via simulations.

Extensions to more than two time-varying treatments

Although we introduce our methods with two longitudinal treatments, our approach can be extended to more than two time-varying treatments in a straightforward fashion. In theory, a fully interacted version of model (2.1) can be formed to include all the main effects of and the interactions thereof. Clinical interests and data sparsity on combinations of treatments may also be used to guide the inclusion of interaction terms into the structural model. Suppose is a collection of causal interaction effects of interest, e.g., b1(t) = a1(t)a2(t), the general joint marginal structural proportional hazards model is where ψ1’s and ψ2’s capture the causal main and interaction effects on the counterfactual hazard function. A consistent estimator of ψ = {ψ11, …, ψ, ψ21, …, ψ2} can be obtained by solving the general form of the estimating equation where Z(A1, …, A, t) is a vector of length W + J representing the time-varying treatment status A(t) and multiplicative terms of the treatment status is evaluated using weighted risk set indicators . The joint treatment weights can be estimated by assuming a specific order in which treatments are initiated and calculating the weights using appropriate history information and , similar as described in Section 3.2. The estimation of the censoring weights Ω(G) also follows the same strategy outlined in Section 3.2 with two longitudinal treatments.

Simulation Study

Simulation design

We carry out simulations to investigate the finite-sample properties of the proposed weight estimators for the marginal structure Cox model parameters. We simulate data compatible with the marginal structural Cox model by generating and relating data adhering to the structural nested accelerated failure time (SNAFT) model (Young et al., 2008). A general representation of a SNAFT model for time-varying treatment a is (Hernán et al., 2005) where is the counterfactual failure time under no treatment. Robins (1992) developed a simulation algorithm to generate data adhering to the SNAFT model under the discrete-time version of the identifying assumptions (A1)–(A3) in Section 3. Young et al. (2008) showed that, under the same identifying assumptions, data adhering to a marginal structural Cox model of the form can be simulated from a SNAFT model with ψaft = ψmsm by adding an additional quantity to the term exp[ψafta(t)]. In particular when has an exponential distribution, the additional quantity is zero, hence the structural nested AFT model simulation algorithm (Robins, 1992) can be used to appropriately simulate data compatible with the marginal structural cox model under complex time-varying data structures. Building on these previous works, we extend the simulation algorithm described in Karim et al. (2017) to generate data from the joint marginal structural Cox model, while allowing for multiple time-varying treatments with discontinuous intervals of treatment eligibility and for both continuous and discrete time-varying confounders.

Throughout we simulate an observational study with n = 1000 patients and two time-varying treatments A1(t) and A2(t). We assume is appropriately summarized by a continuous time-varying confounding variable L1(t) and a binary time-varying confounding variable L2(t). The simulation algorithm includes two steps. In step (1), we consider nonlinear terms for the continuous variable L1(t) and the interaction term A1(t)×L1(t), A2(t)×L1(t), A1(t)×L2(t) and A2(t) × L2(t) in the true treatment decision model. In particular, past treatment status {A1(t), A2(t)} is a predictor of L(t), which then predicts future treatment exposure {A1(t), A2(t)} as well as future failure status Y(t) via (Karim et al., 2017). Therefore, L(t) is a time-dependent confounder affecting both the future treatment choices and counterfactual survival outcomes. The simulation of treatment initiation is placed in the recurrent events framework. Once treatment is initiated at time t, treatment duration following initiation is simulated from a zero-truncated Poisson distribution. In step (1), we generate a longitudinal data set with 100 × 1000 observations (100 aligned measurement time points for each of n = 1000 individuals). In step (2), we randomly discard a proportion of follow-up observations for a randomly selected subset of individuals (Lin et al., 2004); and in the resulting data set, the individuals will have varying number of follow-up measurement time points, which are also irregularly spaced. Web Appendix S2 provides the full pseudo-code for simulating data under the marginal structural Cox model with two time-varying treatments.

Our simulation parameters are chosen so that the simulated data possess similar characteristics to those observed in the COVID-19 data set. The treatments A1 and A2 are simulated to resemble dexamethasone and remdesivir such that: (i) about 20% patients did not take any of the anti-viral and anti-inflammatory medications aimed at treating COVID-19; (ii) among those who were treated, 62% took dexamethasone only, 25% took remdesivir only and 13% took both (either concurrently or with treatment switching); (iii) the number of initiations for both treatments ranges from 0 to 4 with the average medication duration about 5 days. The values of treatment effect parameters ψ1 and ψ2 were set to yield a 6.7% mortality proportion among those who received dexamethasone and a 4.9% mortality proportion in those treated with remedesivir.

Comparison of methods

We conduct two sets of simulations to investigate the finite-sample performance of our proposed joint marginal structural survival model in continuous time (JMSSM-CT). First, we compare how accurate the four weight estimators described in Section 3.3 estimate the structural parameter ψ. Second, we use the best weight estimator, suggested by the first set of simulation, for JMSSM-CT, and compare it with the joint marginal structural model that requires aligned discrete time points (JMSM-DT). To ensure an objective comparison, we use the random forests (Yao et al., 2020) and adapt it into our proposed recurrent events framework to estimate the weights for JMSM-DT. In addition, we implement both JMSSM-CT and JMSM-DT on the “rectangular” simulation data with 100 aligned time points for each individual and on the “ragged” data with irregular observational time points. The performance on the rectangular data will be considered as the benchmark performance, based on which we will assess the relative accuracy of JMSSM-CT and JMSM-DT when estimating the structural parameters with the “ragged” data.

Performance characteristics

To assess the performance of each method, we simulate 250 observational data sets using the above approach, and evaluate the absolute bias, root mean squared error (RMSE) and covarage probability (CP) for estimating the ψ. The CP is evaluated on normality-based confidence intervals with the robust sandwich variance estimator. Figure 2 suggests that the weight estimator (iv) using both the flexible tree-based survival model and kernel function estimator of the treatment initiation intensity yielded the lowest biases in estimating both ψ1 and ψ2. By contrast, the weight estimator (i) via the usual main-effects Cox regression model along with the Nelson-Aalen baseline intensity estimator produced the largest estimation bias. Applying the kernel function smoothing to the Nelson-Aalen estimator led to bias reduction for both the Cox (approach (ii)) and tree-based survival model (approach (iv)) for the treatment process. Flexible modeling of the intensity ratio function has a larger effect in reducing the bias in structural parameter estimates than smoothing the nonparametric baseline intensity estimator. For example, compared to approach (ii), approach (iii) further reduced the mean absolute bias (MAB) in estimating by approximately 67%. Supplementary Table 1 summarizes the MAB, RMSE and CP for the four weight estimators and similarly suggests that approach (iv) led to the smallest MAB and RSME, and provided close to nominal CP with the robust sandwich variance estimator.

Fig. 2.

Biases in the estimates of ψ1 and ψ2 among 250 data replications using four approaches to estimate the weights as described in Section 3.4. Approach (i) uses main-effects Cox regression model and Nelson-Aalen estimator for baseline intensity. Approach (ii) uses kernel function smoothing of the Nelson-Aalen estimator in approach (i). Approach (iii) uses a survival forests model that accommodates time-varying covariates and Nelson-Aalen estimator for baseline intensity. Approach(iv) uses kernel function smoothing of the Nelson-Aalen estimator in approach (iii).

The second set of simulation benchmarks the performance of JMSSM-CT versus JMSM-DT on the data with fully aligned follow-up time points and compare how much each method can recover the benchmark performance in situations where the longitudinal measurements are irregularly spaced. Table 2 displays the MAB, RMSE and CP for each of the two methods under both data settings, and Supplementary Figure 1 visualizes the distributions of biases across 250 data replications. In the rectangular data setting with fully aligned time points, compared to JMSM-DT, JMSSM-CT had similar CP but smaller MAB and RMSE. As the sparsity of longitudinal measurements increased and the time intervals became unevenly spaced, the proposed JMSSM-CT could still recover the benchmark performance; whereas the JMSM-DT had a deteriorating performance (larger MAB and RMSE and lower CP), with larger performance decline under coarser discretization of the follow-up time. Supplementary Table 2 summarizes the distribution of estimated individual time-varying weights from one random replication of the ragged data for JMSSM-CT and JMSM-DT. Overall, JMSSM-CT with the weight estimator (iv) provided the smallest maximum/minimum weight ratio (2.36/0.68) and no extreme or spiky weights.

Table 2.

The joint and interactive effect estimates (log hazard ratio) of COVID-19 treatments and associated 95% confidence intervals (CI), using the COVID-19 dataset drawn from the Epic electronic medical records system at the Mount Sinai Medical Center. The composite outcome of in-hospital death or admission to ICU was used for the general COVID-19 patients. Time to in-hospital death was used for subpopulations of those who had never been admitted into ICU (pre-ICU) and of post-ICU patients. Confidence intervals were estimated via the robust sandwich variance estimators. “×” denotes treatment interaction and “—” indicates that the pairwise interaction was not included in the joint marginal structural survival model.

Treatment classes	ψ^ (95% Confidence Interval)
	Overall	Pre-ICU	Post-ICU
Dexamethasone	−0.20(−0.35, −0.06)	−0.06(−0.71, 0.64)	−0.75(−1.42, −0.08)
Remdesivir	−0.53(−0.75, −0.31)	−0.62(−1.22, −0.02)	−0.54(−1.04, −0.04)
Corticosteroids other than dexamethasone	−0.08(−0.29, 0.19)	−0.13(−0.46, 0.21)	−0.19(−0.27, −0.03)
Anti-inflammatory medications other thancorticosteroids	−0.05(−0.56, 0.47)	−0.28(−1.02, 0.45)	−0.08(−0.89, 0.72)
Remdesivir × Corticosteroids other than dexamethasone	−0.74(−0.95, −0.52)	—	−0.71(−1.38, −0.04)
Dexamethasone × Corticosteroids other than dexamethasone	—	—	−1.13(−1.78, −0.46)

Estimating Causal Effects of Multiple COVID-19 Treatments

We apply the proposed JMSSM-CT with the best weight estimation approach (iv) (Section 3.3) to analyze a comprehensive COVID-19 data set drawn from the Epic electronic medical records system at the Mount Sinai Medical Center, with the goal of studying the comparative effectiveness of multiple COVID-19 treatment strategies. The data set includes 11,286 de-identified unique adult patients (⩾18 years of age) who were diagnosed with COVID-19 and hospitalized within the Mount Sinai Health System between February 25, 2020 to February 26, 2021. A confirmed case of COVID-19 was defined as a positive test result from a real-time reverse-transcriptase PCR-based clinical test carried out on nasopharyngeal swab specimens collected from the patient (Wang et al., 2020). We focus on four treatment classes that are of most clinical interest: (i) remdesivir; (ii) dexamethasone; (iii) anti-inflammatory medications other than corticosteroids; and (iv) corticosteroids other than dexamethasone. Detailed definitions of the four treatment classes are provided in Supplementary Table 3. The observed treatment patterns are visualized in Figure 1; patients could be simultaneously prescribed two or more treatment strategies, or they could switch from one treatment class to another during their hospital stays.

We considered a composite outcome, ICU admission or in-hospital death, whichever occurs first. The outcome may be right censored by hospital discharge or administratively censored on t= February 26, 2021, the date on which the database for the current analysis was locked. For modeling purposes, we assumed that was appropriately summarized by age, sex, race, ethnicity, D-dimer levels (the degradation product of crosslinked fibrin, reflecting ongoing activation of the hemostatic system), serum creatinine levels (a waste product that forms when creatine breaks down, indicating how well kidneys are working), whether the patient used tobacco at the time of admission, whether the patient was admitted to ICU, history of comorbidities represented by a set of binary variables: hypertension, coronary artery disease, cancer, diabetes, asthma and chronic obstructive pulmonary disease, hospital site, and patient oxygen levels (definition provided in Supplementary Table 4). The time-varying confounding variables were ICU admission, D-dimer levels, serum creatbinine level and patient oxygen levels. When fitting the joint marginal structural proportional hazards model (3.12), pairwise treatment interactions were included if there were sufficient data points supporting the joint use of the pair of treatments.

Using the stabilized inverse probability weights to correct for time-varying confounding and censoring, the structural model parameter estimates (log hazard ratio) and the associated 95% confidence intervals are provided in Table 2. Using the parameter estimates, we further computed the counterfactual survival curves under each treatment regimen. Figure 3 presents the counterfactual survival curves for ICU admission or death, whichever occurs first, among patients with COVID-19 infection, under five treatment regimens given upon admission to hospital. Among the four separate treatment classes, remdesivir had significantly better treatment benefits followed by dexamethasone than two alternative treatment classes: anti-inflammatory medications other than corticosteroid and corticosteroids other than dexamethasone. Interestingly, remdesivir and corticosteroids other than dexamethasone had a significant treatment interaction effect suggesting additional survival benefit when they are used in combination. This is demonstrated by the highest counterfactual survival curve under the concomitant use of these two types of medications.

Fig. 3.

Counterfactual survival curves for each of five treatment strategies among the general COVID-19 patients. The composite outcome of ICU admission or death is used.

Discussion

Motivated by inconclusive real-world evidence for the comparative effectiveness of multiple treatment strategies for COVID-19, we have developed a joint marginal structural survival model and novel weighting schemes to address time-varying confounding and censoring in continuous time. There are three main advantages of our proposed method. First, this approach casts the complex time-varying treatment with irregular “start/stop” switches into the process of recurrent events where treatment initiation can be considered under the recurrent event framework with discontinuous intervals of eligibility. This innovative formulation enables us to address complex time-varying confounding by modeling the intensity processes of the filtered counting processes for complex time-varying treatments. Second, the proposed method is able to handle a complex longitudinal dataset on its own terms, without discretizing and artificially aligning measurement times, which would lead to less accurate and efficient treatment effect estimates, as demonstrated by our simulations. Third, modern machine learning techniques designed for censored survival data and smoothing techniques of the baseline intensity function can be used easily for with our weighting method to further improve the treatment effect estimator under conventional parametric formulations. We have also introduced a simulation algorithm that is compatible with the complex data structures of our proposed modeling framework, and demonstrated the accuracy of the proposed method for estimating causal parameters.

Our approach can be extended in the following two directions. First, we considered a joint marginal structural proportional hazards model and a tailored simulation algorithm to generate datasets of complex time-varying structures that are compatible with the proportional hazards model. It may be worthwhile to develop alternative joint marginal structural survival models such as the structural additive hazards model, and assess the robustness of different structural models for estimating counterfactual survival functions under different data generating processes. Second, we have maintained the conditional exchangeability assumption in our work, and therefore developing sensitivity analyses to capture the effects of time-varying unmeasured confounding for our model would be a worthwhile and important contribution.

Table 1.

Comparing the proposed method JMSSM-CT in continuous time with JMSM-DT in discrete time in estimating the treatment effect ψ on the bases of mean absolute bias (MAB), root mean square error (RMSE) and coverage probability (CP) across 250 data replications. In the estimation of the weights, the weight esimator (iv) was used for JMSSM-CT and the random forests adapted into our recurrent events framework (Section 3.2) was used for JMSM-DT. Both methods were implemented on the “rectangular” simulation data with 100 aligned time points for each individual and on the “ragged” data with unaligned time points. With the ragged data, the follow-up time was discretized in the space of 0.5, 1 and 2 days for JMSM-DT.

Data format	Methods	ψ 1			ψ 2
		MAB	RMSE	CP	MAB	RMSE	CP
Rectangular	JMSM-DT	.021	.026	.944	.019	.023	.948
	JMSSM-CT	.015	.020	.948	.014	.018	.948
Ragged	JMSM-DT (2d)	.040	.047	.660	.035	.041	.668
	JMSM-DT (1d)	.033	.041	.732	.029	.035	.738
	JMSM-DT (0.5d)	.027	.034	.801	.024	.030	.804
	JMSSM-CT	.016	.022	.952	.015	.019	.952