SyncTwin: Treatment Effect Estimation with Longitudinal Outcomes
Year: 2021

Abstract:
Most of the medical observational studies estimate the causal treatment effects using electronic health records (EHR), where a patient’s covariates and outcomes are both observed longitudinally. However, previous methods focus only on adjusting for the covariates while neglecting the temporal structure in the outcomes. To bridge the gap, this paper develops a new method, SyncTwin, that learns a patient-specific time-constant representation from the pre-treatment observations. SyncTwin issues counterfactual prediction of a target patient by constructing a synthetic twin that closely matches the target in representation. The reliability of the estimated treatment effect can be assessed by comparing the observed and synthetic pre-treatment outcomes. The medical experts can interpret the estimate by examining the most important contributing individuals to the synthetic twin. In the real-data experiment, SyncTwin successfully reproduced the findings of a randomized controlled clinical trial using observational data, which demonstrates its usability in the complex real-world EHR.

heading: 1 Introduction
text: Estimating the causal individual treatment effect (ITE) using observational data has become increasingly common in the medical literature due to the popularization of electronic health records (EHR). The EHR is a longitudinal collection of records: it contains repeated measurements of a patient’s health condition over irregular time intervals. Although the treatments may also vary over time, many studies consider a point treatment, where the treatment allocation is performed at some observed time and stays fixed during the study [12, 44, 34, 5, 51]. Point treatment is widely applicable to problems involving one-off treatments (e.g. surgical operations) or treatments that do not change frequently (e.g. long term medication for chronic disease). We will refer to the setting above as Longitudinal and Irregularly sampled data with Point treatment setting, or LIP. The LIP setting will be the focus of this work (Figure 1 A).
The LIP setting is different from the conventional settings in causal inference. This is because we observe the pre-treatment outcome yt over time leading to the treatment allocation. These pretreatment outcomes may unveil the inherent temporal structure in the outcome time series (e.g. trend and seasonality), leading to better ITE estimation over time. In contrast, the conventional settings
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
only consider the outcomes after treatment initiation [54, 52, 36, 43, 6]. This is true even for works that deal with dynamic treatment allocations [40, 35, 11].
A naive way to incorporate pre-treatment outcomes into the standard ITE methods is to treat them as additional covariates. This approach is viable but inadequate: the pre-treatment outcomes are arguably much more closely linked to the outcomes after treatment than the covariates — they hence deserve special considerations, e.g. modifying the architecture or loss function to reflect their importance. Without such modification, the ITE methods will not incorporate our prior belief on the importance of the pre-treatment outcomes, which may lead to worse performance.
To bridge this gap, we propose SyncTwin, a novel ITE estimation method tailored for the LIP setting. We assume that the temporal outcomes are generated by individualized latent factors and timevarying latent trends. This assumption is similar to the “factor model” assumption commonly used in Econometrics [1] and allows us to achieve the right balance between the parametric assumption and modeling flexibility. Figure 1 B illustrates the schematics of SyncTwin using an example with two treatment options (treated and control). SyncTwin first uncovers the individualized latent factors using representation learning. For a target individual, SyncTwin selects and weights a few contributors based on their latent factors and a sparsity constraint. It proceeds to construct a synthetic twin whose temporal outcomes are the weighted average of the contributors. Finally, the ITE is estimated as the difference between the target individual and the synthetic twin’s outcomes.
Unlike conventional methods, SyncTwin does not use the learned latent factors to directly predict the outcomes; instead it uses these variables to find the contributors and their importance weights so as to construct the synthetic twin. This approach brings two bonus features; both are important to medical applications. (1) We can calculate an individualized estimation error bound based on the actual and the synthetic outcomes before treatment. In practice, the clinician can accept the recommended treatment when the error bound is below a threshold and resort to expert knowledge otherwise. (2) We can identify the most important contributors to an estimate as the ones who receive the highest weights. The clinician can further examine these contributors (e.g. whether they are indeed similar to the target individual) to validate and interpret the estimate.
Contributions. (1) We develop SyncTwin to leverage the temporal structure in the outcomes for better ITE estimation. (2) We provide theoretical justifications for each step involved in SyncTwin and prove an individualized error bound that can be used to identify untrustworthy estimates. (3) In addition to extensive simulations, we conduct an observational study using real EHR data and successfully reproduced the findings of a randomized controlled clinical trial with SyncTwin.

heading: 2 Problem setting
text: We consider an observational study with N individuals indexed by i ∈ [N ] = {1, . . . , N}. Let ai ∈ {0, 1} be the treatment indicator with ai = 1 if i received the treatment at some time and ai = 0
otherwise1. We realign the time steps such that all treatments were initiated at time t = 0. Let I1 = {i ∈ [N ] | ai = 1} and I0 = {i ∈ [N ] | ai = 0} be the set of the treated and the control individuals respectively. Denote N0 = |I0| and N1 = |I1| as the sizes of the groups. Let Xi = [xis]s∈[Si] be the temporal covariates consisting of a sequence of observations xis ∈ RD. Let Mi = [mis]s∈[Si] be the sequence of binary masking vectors mis ∈ {0, 1}D, where [mis]d = 1 indicates the dth element of xis is measured and [mis]d = 0 otherwise. The entire sequence Xi contains Si ∈ N observations taken before treatment at times Ti = [tis]s∈[Si], tis < 0. The maximum sequence length S = maxi(Si).
The outcome yit ∈ R is observed at time t ∈ T − ∪ T +. Let T − = {−M, . . . ,−1} and T + = {0, . . . ,H − 1} be the observation times before and after treatment allocation. We arrange the outcomes after treatment into a H-dimensional vector denoted as yi = [yit]t∈T + ∈ RH . Similarly define pre-treatment outcome vector y−i = [yit]t∈T − ∈ RM . Using the potential outcome framework [41], let yit(ai) ∈ R denote the potential outcome at time t in a world where i received the treatment as indicated by ai. Let yi(1) = [yit(1)]t∈T + ∈ RH , and y−i (1) = [yit(1)]t∈T − ∈ RM . Similarly, let yi(0) = [yit(0)]t∈T + and y − i (0) = [yit(0)]t∈T − . The individual treatment effect (ITE) is defined as τi = yi(1)− yi(0) ∈ RH .

heading: 3 ITE estimation via SyncTwin
text: SyncTwin estimates the ITE by predicting both potential outcomes yi(1) and yi(0) and take the difference. Note that for an individual in the dataset i ∈ [N ], it is only necessary to predict the counterfactual outcome yi(1− ai) because the factual outcome yi(ai) = yi is observed (under the consistency assumption to be discussed later). Without loss of generality, we will use estimating yi(0) as an example in the following sections. One can estimate yi(1) using the same method.

heading: 3.1 Assumptions
text: SyncTwin relies on three assumptions. (1) Stable Unit Treatment Value Assumption [41]: yit(ai) = yit, ∀i ∈ [N ], t ∈ T − ∪ T +. (2) No anticipation, also known as causal systems [4]: yit = yit(1) = yit(0), ∀t ∈ T −, i ∈ [N ]. (3) Data generating process (DGP). The DGP assumption involves two parts. First, the causal directed acyclic graph (DAG) is specified in Figure 2 [37], which involves a latent factor ci ∈ RK . Secondly, we assume the potential outcomes have a parametric form:
yit(0) = q > t ci + ξit, ∀t ∈ T − ∪ T +, (1)
where qt ∈ RK , K < min(M,H) is a weight vector and ξit is the white noise. Equation 1 is commonly referred to as a latent factor model with ci as the individualized latent factor and qt as the time-varying latent trend [10]. Without loss of generality, we require the weight vectors ||qt|| = 1, ∀t ∈ T − ∪ T + [50].
Discussion. SyncTwin differs from the nonparametric ITE estimation methods [27] because it assumes the parametric form in Equation 1. In essence, it means that we can fully separate the time effect from the individual effect. We would like to highlight two aspects of this assumption. First, the latent factor ci does not change over time. It provides a constant link between the pre-treatment outcomes and the post-treatment outcomes. Secondly, the time-varying latent trend qt is common to all individuals. Together, they ensure that any linear combination of different individuals’ outcomes will automatically preserve qt, i.e. ∑ i biyit(0) = q > t ∑ i bici + ∑ i biξit, for any weights bi ∈ R, ∀i ∈ [N ]. For this reason, SyncTwin
bypasses qt and directly finds the weights bi’s to issue counterfactual predictions. As we will show later, doing so allows us to derive various theoretical results to inform the design of SyncTwin . It also leads to a checking procedure to validate estimation quality and improvements in interpretability.
1SyncTwin can model more than two treatment options, but we focus on binary treatments for illustration.
The DAG in Figure 2 is different from the one commonly used in the static setting [41]. We assume the latent factor ci captures the physiological factors that impact both the outcome and the covariates. Here, the covariates xis is a “confounder” in a general sense because it opens up a backdoor path from treatment to outcome, i.e. ai → xis → ci → yis [37]. Based on the backdoor criterion [37], adjusting for the covariates xis is sufficient to identify the treatment effect from observational data.
We compare our assumptions with those used in the related works in Appendix A.2. In particular, we show that our assumption is weaker than the one used by Synthetic Control, a widely-applied method in Econometrics [2]. In Appendix A.5, we further demonstrate the plausibility of our assumptions in real-world scenarios. In the simulation study (Section 5.1) we show experimentally that SyncTwin performs well even when the data is not generated from the assumed DGP exactly but instead from a set of differential equations. Finally, and perhaps most importantly, we show that based on our assumptions, SyncTwin is able to reproduce the findings of a large-scale randomized controlled trial using observational data (Section 5.2).

heading: 3.2 Learning to represent temporal covariates
text: The latent factor ci plays an important role as it affects the covariates xis and the outcomes yit. SyncTwin uses deep neural networks to learn the representation c̃i as a proxy for ci.
SyncTwin is agnostic to the exact choice of architecture as long as the network translates the irregularly sampled temporal covariates into a fixed-sized representation vector. For this reason, we use the well-proven sequence-to-sequence architecture (Seq2Seq) [47] with a standard attentive encoder [9] and a LSTM decoder [26]. The learned representation c̃i = fe(Xi,Mi, Ti; θe), where fe is the encoder with trainable weights θe. The reconstructed covariates X̂i = fd(c̃i, Ti; θd), where fd is the decoder with trainable weights θd. To ensure the learned representation c̃i is a linear predictor of the potential outcome yi(0) (Equation 1), we introduce a trainable parameter Q̃ ∈ RH×K and define ỹi(0) := Q̃ · c̃i. Note that using a nonlinear function to map from c̃i to ỹi(0) will be inconsistent with the DGP and will not uncover the latent factor ci as desired.
We train the networks end-to-end by optimizing the loss function Ltr = λrLr + λpLs, where λr and λp are hyper-parameters2 balancing the supervised loss Ls and the reconstruction loss Lr:
Ls(D0) = ∑ i∈D0 ||ỹi(0)− yi(0)||2; Lr(D0,D1) = ∑ i∈D0∪D1 ||(X̃i −Xi) Mi||2, (2)
where D0 ⊆ I0 and D1 ⊆ I1 are training data, mis is the masking vector, represents element-wise product and || · || is the L2 norm. In Proposition 1, we show that minimizing the supervised loss Ls will reduce the error bound on the learned representations, making them closer to the true latent factor ci. The proof and the motivations for Lr are shown in A.1.1. Proposition 1 (Error bound on the learned representations). Under the assumptions in Section 3.1, the total error on the learned representations for the control is bounded by:∑
j∈I0 ‖cj − c̃j‖ ≤ βLs + ∑ j∈I0 ‖ξj‖, (3)
where Ls is the supervised loss, β is a constant depending on qt and Q̃, and ξj is the white noise (Equation 1 and 2).

heading: 3.3 Synthesizing the twin
text: At this point, one may be temped to predict the counterfactual outcome yi(0) of a treated individual i ∈ I1 with the output of the neural network ỹi(0). However, ỹi(0) is issued by a black-box neural network, which is not easily interpretable. As a remedy, SyncTwin explicitly constructs a synthetic twin who matches the target in representations. As we will show, this approach is able to control estimation bias and is more interpretable.
The synthetic twin of a target individual is defined by a set of weights bi := [bij ]j∈I0 ∈ RN0 , each associated with a contributor in the control group j ∈ I0 (or treatment group when predicting ŷi(1)).
2The hyperparamter sensitivity is studied in A.13. λr and λp do not significantly impact the performance.
SyncTwin solves the following optimization problem to find the weights:
bi = argmin b̃i ‖c̃i − ∑ j∈I0 b̃ij c̃j‖2 s.t. b̃ij ≥ 0, ∀j ∈ I0 and ∑ j∈I0 b̃ij = 1, (4)
where c̃j is the representation learned by the encoder. Denote the loss function in Equation 4 as Lm. We will call Lm the matching loss because it indicates how well the synthetic twin matches the target individual in representations. Let dci denote the optimal value of Lm and let ĉi be the synthetic representation given by the solution bi:
dci := ‖c̃i − ∑ j∈I0 bij c̃j‖2; ĉi := ∑ j∈I0 bij c̃j (5)
The optimization procedure is detailed in Appendix A.8. SyncTwin predicts the potential outcome yit(0) for a target individual i using the same weights bi: ∀t ∈ T − ∪ T +,
ŷit(0) = ∑ j∈I0 bijyjt(0) = ∑ j∈I0 bijyjt, (6)
where the last equality follows from the consistency assumption (Section 3.1). Denote ŷi(0) = [ŷit(0)]t∈T + as the predicted potential outcomes after treatment initiation.
The estimator in Equation 6 follows directly from the DGP assumption (Equation 1). To see this, remember that yit(0) is linear in the latent factor ci. Hence, if we find a set of weights b ∗ ij to match the latent factor, i.e. ci = ∑ j b ∗ ijcj , the same set of weights will also match the outcome
yit(0) ≈ ∑ j b ∗ ijyjt(0) up to random noise. In practice, since we do not observe ci, we have to perform matching on the learned representation c̃i as in Equation 4.
Therefore we can see that the quality of the estimator ŷit(0) depends on two things: (1) whether the learned representation c̃ is close to the true latent factor c and (2) whether good matching weights bi can be found on the learned representation. SyncTwin attempts to fulfill the first condition by representation learning in Section 3.2. The second condition is facilitated by solving the optimization problem in Equation 4. Proposition 2 formalizes this intuition and shows that when a perfectlymatching twin is found, the estimation error only depends on the quality of the learned representations (proved in A.1.1). Proposition 2 (Bias bound on counterfactual prediction). Suppose that dci = 0 for some i ∈ I1 (Equation 1 and 5), the absolute value of the counterfactual prediction bias for i is bounded by:
|E[ŷi(0)− yi(0)]| ≤ |T +| ( ‖ci − c̃i‖+ ∑ j∈I0 ‖cj − c̃j‖ ) .
Notes on Interpretability. In addition to reducing bias by learning and matching representations, SyncTwin is also more interpretable. In particular, the weights bi can be interpreted as the “contribution” or “importance” of a contributor j to the target i due to the optimization constraints. We can create a shortlist of the most important contributors based on bi. A domain expert can understand the rationale behind the estimate by examining the shortlist. For instance, they can check whether the important contributors share similar disease progression patterns as the target individual.
This procedure corresponds to the notion of data point interpretability or example-based interpretability [45], where one explains the prediction by presenting the most relevant data points to the users. Note that the learned representations c̃i are internal to the algorithm; we do not intend to show these representations to the user or to make interpretations on them.

heading: 3.4 Calculating individualized error bound with pre-treatment outcomes
text: Since Equation 6 applies to all time points before and after treatment, let ŷ−i (0) = [ŷit(0)]t∈T − be the predicted pre-treatment potential outcome vector. We define dyi as the estimation error in the pre-treatment period:
dyi = ‖ŷ − i (0)− y − i (0)‖1 = ‖ŷ − i (0)− y − i ‖1, (7)
where || · ||1 is the vector `1-norm. The second equality allows dyi to be evaluated for all individuals. It holds because of the no anticipation assumption (Section 3.1). Achieving a small error dyi implies
that the synthetic twin ĉi matches well with the true latent factor ci due to the linearity between y−i (0) and ci (Equation 1). Since ci is assumed to be constant over time, the pre-treatment error can be used to assess the post-treatment error. We formalize this intuition in Proposition 3 and use dyi to control the error in the post-treatment period (proved in A.1.1). Proposition 3 (Error control under no hidden confounders). Given any target error threshold δ > 0, define the acceptance group of treated individuals as
Aδ = { i ∈ I1|dyi ≤ δ|T −|/|T +| } .
Under the assumptions in Section 3.1, the post-treatment estimation error |E[ŷi(0)]− E[yi(0)]| ≤ δ, ∀i ∈ Aδ .
Proposition 3 shows that we can control the estimation error to be below a certain threshold δ by rejecting the estimate if its error dyi during the pre-treatment period is larger than δ|T −|/|T +|. Alternatively, we can rank the estimation trustworthiness for the individuals based on dyi . This is helpful when the user is willing to accept a percentage of estimations which are deemed most trustworthy.

heading: 3.5 Training, validation and inference
text: We perform model training, validation and inference (testing) on three disjoint datasets. In summary, we train the encoder and decoder on the training data using the loss functions described in Section 3.2. The validation data is then used to validate and tune the hyper-parameters of the encoder and decoder. Then, we solve the optimization problem in Section 3.3 to find the weights bi for the target individuals in the testing data. Finally, we compute the potential outcome estimator (Section 3.3) and the individualized error bound (Section 3.4). The pseudocode is described in A.7.

heading: 4 Related work
text: Synthetic Control. Similar to SyncTwin, Synthetic control (SC) [1] and its extensions [8, 7] estimate ITE based on synthetic outcomes. However, SC needs to flatten the temporal covariates [xis]s∈[Si] into a fixed-sized (high-dimensional) vector xi and use it to construct the twin. As a result, SC does not allow the covariates to be variable-length or sampled at irregular frequencies (Table 1); otherwise xi’s dimensionality will vary across individuals. This severely limits SC’s usability in medical observational studies because such sampling irregularities are common in EHR. Moreover, SC assumes yit(0) = q > t xi + ξit, i.e. the flattened covariates xi linearly predicts yit(0), which is a special case of SyncTwin’s assumption (to be discussed later in Equation 1) and seems unlikely in many medical applications (e.g. body mass index does not linearly relate to blood pressure).
Static nonparametric ITE estimation. Most ITE estimation methods in the literature consider the static setting and do not pay special attention to the outcomes before treatment y−i — these outcomes are treated as standard covariates Xi, if at all (Table 1). However, the pre-treatment outcomes y−i may encode certain temporal structure that will persist even on the outcomes yi after treatment (e.g. trend and seasonality). Explicitly unveiling such temporal structure will be very beneficial to accurately estimating and validating the ITE after treatment. Moreover, these methods make no parametric assumption on the data generating model [27] and they learn to predict the two potential outcomes
using deep neural networks [43, 54, 22]. Despite the recent progress in interpretable machine learning, these methods are still in general less interpretable than the parametric models such as SyncTwin whose coefficients can be directly understood by an expert.
ITE estimation for dynamic treatments. Several works consider ITE estimation for dynamic treatments (DT) [40]. In the DT setting (Table 1), multiple treatment decisions are made over time for each individual. ITE estimation is much more challenging for DT because a prior treatment at−1 may influence the temporal covariates xt, which may in turn confound a latter treatment at. This problem is known as temporal confounding. Marginal structural model [40] and related works [42, 35, 11] have been proposed to overcome temporal confounding. In comparison, our setting considers a point treatment and therefore do not suffer from temporal confounding. As a result, although the works in the DT setting are applicable, they might be over-complicated for our problem. In addition, our previous discussion about the inadequacy in leveraging pre-treatment outcomes also applies to DT methods.
Works with similar terminology. Several unrelated works in the literature use terms such as “twin”, which might cause confusion. We discuss these works in Appendix A.4.

heading: 5 Experiments

heading: 5.1 Simulation Study
text: Data generation.3 In this simulation study, we evaluate SyncTwin on the task of estimating the LDL cholesterol-lowering effect of statins, a common drug prescribed to hypercholesterolaemic patients. Following our convention, the individuals are enrolled at t = 0, the covariates are observed in T = [−S, 0), where S ∈ {15, 25, 45}, and the ITE is to be estimated in the period T + = [0, 4]. We start by generating a covariate kint for from a mixture distribution:
kinit = δifit + (1− δi)git; δi iid∼ Bern(pi) (8)
where Bern(p) is the Bernoulli distribution with success probability p and fit, git are drawn from two different distributions (specified in A.10). To introduce confounding bias, we vary pi for the treated and the control: pi = p0, ∀i ∈ I0 and pi = 1, ∀i ∈ I1. The constant p0 controls the degree of confounding bias (smaller p0, larger bias). We simulate the outcome yt using the widely adopted Pharmacology model in the literature [16, 53, 32]. It is obtained by solving the differential equation 9 and adding independent white noise ∼ N(0, 0.1).
ṗt = k in t − k · pt; ḋt = at − h · dt; ẏt = k · pt − kdt(dt + d50)−1 · yt. (9)
where ẏt is the derivative of yt and at is the indicator of statins treatment. The covariates includes xt = {kint , yt, pt}. The interpretation of all constants involved are presented in Appendix A.10. Finally, we introduce irregular sampling by creating masks mit∼Bern(m), where probability m ∈ {0.3, 0.5, 0.7, 1}. It is worth highlighting that the simulation data are not generated from SyncTwin’s assumed DGP (1) but from domain-specific ODEs (9).
Benchmarks. From the Synthetic Control literature, we considered the original Synthetic Control method (SC) [2], Robust Synthetic Control (RSC) [7] and MC-NNM [8]. From the dynamic treatment literature, we compared against Counterfactual Recurrent Network (CRN) [11] and Recurrent Marginal Structural Network (RMSN) [35]. Note that RMSN belongs to the family of Marginal Structural Models (MSMs) and we use it as a representative of MSMs. We included a modified CFRNet, which was originally developed for the static setting [43]. We replaced its fully-connected encoder with the encoder used by SyncTwin to to allow CFRNet to model temporal covariates (Section 3.2). We also included a benchmark adapted from the counterfactual Gaussian Process (CGP) [42] in order to adjust for the patient level covariates. We use One-nearest Neighbour Matching (1NN) [46] as a baseline. The implementation details of all benchmarks are available in Appendix A.9. We compared with two ablated versions of SyncTwin. SyncTwin-Lr is trained only with reconstruction loss and SyncTwin-Ls only with supervised loss.
3The implementation of SyncTwin and the experiment code are available at https:// github.com/ZhaozhiQIAN/SyncTwin-NeurIPS-2021 or https://github.com/orgs/ vanderschaarlab/repositories
Evaluation metric. In many medical studies, the outcome of interest is measured with heavytailed non-Gaussian noise [30]. In these settings, the mean absolute error (MAE) is preferred over the mean squared error (MSE) as an evaluation metric due to its robustness to the outliers [38]. For this reason, we will evaluate the mean absolute error (MAE) on counterfactual prediction: 1 N1 ∑N1 i=1 ||yi(0) − ŷi(0)||1. Note that MAE is also directly comparable to the error bound in Proposition 3 (i.e. the L1 norm).
Main results. Table 2 presents the results for various levels of confounding bias p0. Additional results for different sequence length S and sampling irregularity m are shown in the tables in Appendix A.11. SyncTwin achieves the best or equally-best performance in all cases despite the assumed DGP does not exactly hold. The full SyncTwin with both loss functions also consistently outperforms the versions trained only with Lr or Ls. In comparison, SC, RSC and MC-NNM underperform because their assumption that the flattened covariates xi linearly predict the outcome is violated (Section 4).
Error bound. SyncTwin allows the user to reject untrustworthy estimates based on dyi . In Figure 3 (A) below, we show the scenarios of using various rejection thresholds on dyi . As expected from Proposition 3, when the user increases the threshold, more estimates will be accepted but they will have a higher MAE. The fact that the two curves in the figure share the similar increasing trend suggests dyi is a good indicator of individual error: rejection based on a non-informative criterion will lead to a flat line of MAE. In practice, the user can use Figure 3 (A) to calibrate the threshold in order to achieve balance between the accuracy and the workload. For instance, if the user would like SyncTwin to accept 75% of the estimate, he or she should set threshold δ = 0.12 according to the pre-treatment error dyi and expect a MAE of around 0.105 for the post-treatment outcome.
Sensitivity to omitted covariates. In Figure 3 (B), we show the situation when some covariates are omitted from the analysis, thus making it harder to adjust for the confounding bias. In particular,
we set xt = {kint , yt, pt}, xt = {yt, pt} or xt = {yt} (omitting 0,1,2 covariates). We observe that SyncTwin is still the best-performing model when some covariates are omitted.

heading: 5.2 Experiment on real data
text: Few existing works in the literature validate the method by conducting a real-world observational study to replicate the findings of a randomized controlled trial (RCT), which is the gold standard of causal inference. Most existing works rely on synthetic or semi-synthetic data for validation, which may significantly under-represent the complexity of the real-world data.
To validate the usability of SyncTwin in real-world medical problems, we conduct an observational study to emulate a large-scale RCT—Heart Protection Study (HPS) [21, 20]. Although most existing methods have been validated in synthetic or semi-synthetic experiments [43, 11], few of them have been demonstrated to successfully reproduce the findings of a RCT from observational data. Since RCT is the gold standard of treatment effect estimation, the ability to reproduce a RCT provides strong evidence about the method’s usability in medical research and facilitates the method’s adoption in the medical community.
The HPS was conducted to investigate the treatment effect of statins, a drug commonly used to lower the LDL Cholesterol (LDL). It reported an a change of -1.26 mmol/L (SD=0.06) in LDL for participants randomised to statins versus placebo at the end of the first year after treatment [39]. The study enrolled 20,536 individuals and lasted for eight years, making it one of the largest studies to investigate the treatment effect of statins.
Data Source. We used medical records from English National Health Service general practices that contributed anonymised primary care electronic health records to the Clinical Practice Research Datalink (CPRD), covering approximately 6.9 percent of the UK population [25]. CPRD was linked to secondary care admissions from Hospital Episode Statistics, and national mortality records from the Office for National Statistics. We defined treatment initiation as the date of first CPRD prescription and the outcome of interest was the measured LDL. Known risk factors for LDL were selected as temporal covariates measured before treatment initiation: HDL Cholesterol, Systolic Blood Pressure, Diastolic Blood Pressure, Body Mass Index, Pulse, Creatinine, Triglycerides and smoking status. Our analysis is based on a subset of 125,784 individuals who met the enrollment criterion of HPS (Appendix A.15). They were split into three equally-sized subsets for training, validation and testing, each with 17,371 treated and 24,557 controls.
Main findings. SyncTwin estimates the average treatment effect as -1.25 mmol/L (SD 0.01) by averaging the estimated ITE on testing data ∑ i∈Dte τ̂it/|Dte|. The estimate is very close to the result reported in HPS, -1.26 mmol/L (SD=0.06), after taking into account the standard errors. In contrast, CRN and RMSN estimate the effect to be -0.72 mmol/L (SD 0.01) and -0.83 mmol/L (SD 0.01) respectively, which are significantly smaller than the trial’s finding. Other benchmark methods (e.g. SC, RSC and MC-NNM) either cannot handle irregularly-measured covariates or do not scale to the size of the dataset. Our result suggests SyncTwin is able to overcome the confounding bias in the complex real-world datasets. Note that SyncTwin is also able to estimate ITE and the effect over different time horizons. However, since these results are not reported in HPS, we are not able to quantitatively evaluate them.
Qualitative evaluation of ITE over time. For each individual, we can visualize the outcomes before treatment and compare them with the synthetic twin to sense-check the estimate. The individual shown in Figure 4 (top) has a sensible ITE estimate because the synthetic twin matches its pretreatment outcomes closely. In addition to visualization, we can calculate the individualized error bound based on dyi (Equation 7). From Figure 4 (bottom left) we can see in most cases d y i is small with a median of 0.24 mmol/L (compared to the average between any two randomly chosen individuals 0.76 mmol/L). This means if the expert can only tolerate an error of 0.24 mmol/L on ITE estimation, half of the estimates (those with dyi ≤ 0.24 mmol/L) can be accepted (Section 3). As shown in Figure 4 (bottom right) on average only 15 (out of 24,557) individuals contribute to the synthetic twin (bij > 0.01). The medical experts can check these few contributors to interpret the estimate.

heading: 6 Conclusion and future work
text: We present SyncTwin, an ITE estimation method tailored for temporal outcomes with point treatment. SyncTwin achieves interpretability and strong performance in both simulated and real data experiments. In future works, we plan to extend SyncTwin to the dynamic treatment setting and to model the outcomes in continuous time. Developing additional assumptions to guarantee the robustness to certain types of unobserved confounders is also an interesting avenue for future research.
Acknowledgments and Disclosure of Funding
We thank anonymous reviewers as well as members of the vanderschaar-lab for many insightful comments and suggestions. This work is supported by the US Office of Naval Research (ONR), the National Science Foundation (NSF Grant number:1722516), the Alan Turing Institute (under the EPSRC grant EP/N510129/1), and GlaxoSmithKline (GSK).

Review summaries:
This paper proposes a method, SyncTwin, to estimate treatment effects in Longitudinal and Irregularly sampled data with Point treatment (LIP) setting. The idea is to create a synthetic twin from controls that closely matches a treated target in representation. By using the pretreatment time samples of the twin to estimate temporal trend, SyncTwin can predict the outcome of the target when no treatment is administered. The method is applied to simulated and real data with lower prediction error shown. 

The paper introduces a latent variable model (SyncTwin) for estimating time-varying individualized (point-) treatment effects (ITE) given temporal irregularly sampled pre-treatment covariates and outcomes. Theoretical results justifying the proposed approach include upper bound errors on the learned representations and counterfactual predictions. Experimental results on synthetic and real-world datasets show SyncTwin outperforms baselines per mean absolute error (MAE) on ITE and recovering RCT findings, respectively.

The paper proposes a representation learning approach to model counterfactuals in longitudinal data with point interventions over time. Their approach, SyncTwin is thus suited to situations where an intervention is made on a temporal dataset at an instance of time or as a continuous exposure by treating all points post the intervention as the potential outcomes under treatment.   They compare the performance of their model against several counterfactual estimation baselines originally proposed to model static observational data with interventions.

The paper proposes an individual treatment effect estimation method called SyncTwin, which learns synthetic twins of target patients via the use of representation learning for counterfactual predictions. The proposed method is used on longitudinal data such as electronic health records. The paper is concerned with the setting of longitudinal and irregularly sampled data with point treatment. The paper describes the proposed method, offers some theoretical guarantees, and demonstrates the utility of the method on both synthetic and real-world data.