The study population included a retrospective cohort of patients who were at high risk for stroke and had been treated and followed up at Ramathibodi Hospital, Thailand, between 2010 and 2020. Hospital records and the International Classification of Diseases, 10th Revision (ICD-10) classification system were used to identify patients. Patients were eligible if they were aged >18 years and had one or more of the following conditions: hypertension (HT; ICD-10 code I10-I16), diabetes mellitus (DM; ICD-10 code E08-E13), dyslipidemia (DLP; ICD-10 code E78), and atrial fibrillation (AF; ICD-10 code I48). Patients were excluded if they had a stroke on their first visit or only had one visit during the study period. The main outcome measured in the study was the occurrence of ischemic or hemorrhagic stroke, which was identified using the ICD-10 codes I63 and I61, respectively.

Patients were followed up from their index date (i.e., the date they were identified as high-risk patients) until they progressed to stroke, were lost to follow-up, or were stroke-free at the end of the study (December 31, 2020). Patients who were lost to follow-up or stroke-free at the end of the study period were censored on their last visit date or at the end of the study. A causal diagram was constructed (Figure 1), and potential predictors of stroke were collected, including age, sex, BMI, chronic kidney disease (CKD), AF, HT, DM, and DLP. HT, AF, and DM were considered as mediators, whereas the remaining variables were covariates in the models. A software library called DoWhy, now incorporated into PyWhy (Python Software Foundation), was used to construct models for stratification, IPW, DRE, and DML [24]. Parameters of all estimators were set by default in the DoWhy package. The number of strata in the stratification method was automatically determined [25]. The weighting scheme in IPW was set to default inverse propensity score. For DRE, the regression and propensity models were specified as lasso and logistic regression, respectively. For DML, linear and nonlinear cross-fitted models were applied to the outcome model (lasso and Extreme Gradient Boosting [XGBoost]), propensity model (logistic regression and XGBoost), and final model (linear regression and lasso). Estimands of each risk pathway were defined by PyWhy from the input causal graph. Graphical causal model–based inferences from the DoWhy library were used for medication analysis to quantify the causal effects of direct and indirect pathways, termed natural direct effect (NDE) and natural indirect effect (NIE), respectively [4,26]. NDE (Y1,M(0)x−Y0,M(0)x) refers to the change in the outcome of an individual when they are exposed to a specific treatment Y1, compared to another treatment Y0, while keeping the mediator variable constant at the baseline value or reference treatment M(0). In contrast, NIE (Y1,M(1)x- Y1,M(0)x) refers to the difference between the counterfactual outcome value when treatment Y1 is fixed and the mediator assumes a certain value at a particular treatment M(1) and the counterfactual outcome value when the mediator assumes the same value at the baseline M(0) [27].

The Dragonnet NN was used to estimate PO and propensity scores. The architecture of Dragonnet was based on previous work (Figure 2) [18]. It employs a deep net to create a representation layer ɸ(X) ∈ ℝᴾ, which is used to forecast outcomes for both the treatment Ŷ(1) and control groups Ŷ(0). It utilizes 2 hidden layers for each outcome model while a basic fully connected layer with a sigmoid function is used for the propensity score (𝜺). CATE was estimated by subtracting treatment (risk) and control PO for each risk factor (Y1x−Y0x∣Z) and risk ratios were estimated by division of PO (Y1xY0x| Z); Y₁ is the PO for the risk group, Y₀ is the PO for the control group, x is an interested factor, and Z are other covariates.

To accurately estimate the ITE, it is mandatory for the conditional independence assumption to hold, especially considering the unequal distribution of covariates between factual and counterfactual outcomes of the treatment and control groups, commonly known as covariate shift. To address this challenge, we employed a nested method of weighted split-conformal quantile regression (CQR) to estimate the ITE [20,23] by incorporating antiplatelet medications as a treatment for stroke prevention. POs were estimated using quantile loss setting α at .05. The dataset was split evenly into training and evaluation sets; Multimedia Appendix 1 shows the entire algorithm. All risk factors and covariates were similar between models, considering antiplatelet medication as a treatment and stratified by risk factor (Yantiplatelets=1x- Yantiplatelets=0x| Z), with x representing the risk factors of interest (i.e., HT, DM, and DLP) and Z representing other covariates. AF was not included as a stratum for the estimation of ITE in this example since it is not an indication for the prescription of antiplatelet therapy, but it remained a covariate.

The data were anonymized to ensure confidentiality and privacy protection. This study was approved by the Human Research Ethics Committee, Faculty of Medicine Ramathibodi Hospital, Mahidol University (COA. MURA2021/255). The committee waived the need to obtain consent for the collection, analysis, and publication of the retrospectively obtained and anonymized data for this noninterventional study.

We estimated the causal influences of risk factors associated with stroke outcomes using multiple approaches that included SPS, IPW, DRE, SCM, and mediation analysis, in addition to DML and Dragonnet NNs. Our findings indicate strong positive causal effects associated with AF and HT on stroke development, with DM exerting a weaker effect. DLP, in contrast, had little effect. Furthermore, our analysis suggests that patients with both DM and HT not currently in receipt of antiplatelet treatments would be the most likely beneficiaries of antiplatelet therapy based on the mean ITEs.

The results from the different estimators generally demonstrated consistency, although there were slight variations in specific point estimates and confidence intervals varied slightly. The estimated causal effect derived from various methods using real-world observational data is comparable with standard cohort epidemiological studies using more traditional logistic regression approaches [28,29].

SPS is a widely used method that minimizes confounding bias by adjusting baseline covariates and confounding factors and estimating treatment effects by stratum. However, SPS is sensitive to the number of strata and features that affect both treatment and outcome (confounding factors), which can lead to bias in the causal effect estimate [30-33]. In addition, some strata may be sparsely populated, making the ATE hard to define and more prone to bias [34]. Rosenbaum and Rubin [9] originally proposed dividing the strata into 5 levels and then subsequently automatically splitting the strata until the balance in the numbers of treated and control observations was achieved [25].

IPW attempts to reduce confounding of the ATE by weighting the sample with the inverse propensity score and by balancing the distribution of the covariates between the treated and untreated groups [35], thereby avoiding the problem of data sparsity that may be present in SPS, particularly with small sample sizes. However, there is a reliance on the assumption that the propensity score model correctly captures all confounding factors, which, if incorrect, may bias the ATE. Additionally, IPW is more sensitive to the model and variable selection for estimating the propensity scores, with small differences in estimated propensity scores potentially leading to large differences in estimated causal effects [36]. Finally, IPW may imprecisely estimate treatment effects if a sample size is small, leading to a propensity score close to 0 or 1 [36,37].

DRE combines propensity score and outcome regression models [38], which can lead to improvements in the robustness of model specification by allowing one of the two treatment and outcome models to be miss-specified but still provide a consistent estimation [39]. The challenge is to validly model either the propensity score or the outcome model; it may be tempting to use modern machine learning approaches or nonparametric models in DRE, but this may lead to bias if the functions are too complex, leading to overfitting [40,41]. DML was developed to address the bias from regularization and overfitting in estimating the parameter of interest, which arises when naively inserting machine learning estimators into the estimation equation. This approach consists of two critical components: (1) the use of Neyman-orthogonal moments or scores to estimate the parameters and (2) the application of cross-fitting, which provides an efficient form of data-splitting. By using both elements, DML minimizes the impact of regularization bias and overfitting on parameter estimation; this also extends to nonparametric models [14].

Applying POMs (eg, SPS, IPW, DRE) relies heavily on the assumption that the treatment assignment is independent of the PO given the observed covariates, which is known as “unconfoundedness” or the conditional independence assumption. If this assumption does not hold, the estimated causal effect will be biased. In contrast, SCMs facilitate the modeling of complex relationships between multiple causes and effects in the presence of latent or unobserved variables [4,42]. In addition, SCMs can be considered as counterfactual predictions of interventions, which can be useful in applications such as causal inference in experimental or observational studies [43-46]. However, SCMs are limited by the assumption of independence between variables and may require conceptualized causal relationship mechanisms.

The benefit of using NNs to estimate causal effects is their flexibility and power to handle high-dimensional and complex data. Shalit et al [17] introduced TARNet by sharing information between the PO of treatment and control groups, which is different from the previous model that separated the training data. More recently, Dragonnet was developed by combining propensity scores with targeted regularization, resulting in more accurate inference [18]. Dragonnet is considered more robust with very low or high propensity scores but has several limitations including sensitivity to choice of architecture and hyperparameters, dealing with only a single set of features at a time, and difficulty of interpretation [18]. Despite some limitations, Dragonnet’s benefits surpass these drawbacks, making it an attractive approach for estimating causal effects in complex real-world data.

A critical aspect of causal inference, particularly in estimating CATE, involves certain assumptions, notably ignorability and positivity. Strong ignorability necessitates the observation and adjustment for all confounding variables that influence both the treatment and the outcome, while positivity ensures that every patient has a nonzero probability of receiving each treatment. In our study, we believe these assumptions are reasonably satisfied. We included a comprehensive set of covariates, such as age, sex, BMI, chronic kidney disease, and relevant comorbidities (HT, DM, DLP, and AF), which are well-documented factors influencing stroke risk and treatment decisions. However, we acknowledge that there might be unmeasured confounders not captured in our dataset. Regarding the decision on antiplatelet drug administration, we utilized detailed patient records from Ramathibodi Hospital, ensuring a thorough assessment of factors influencing treatment. Nonetheless, we recognize the potential for residual confounding and the inherent limitations of observational data. Future studies could benefit from incorporating more granular clinical data and leveraging advanced causal discovery methods to further validate these assumptions.

Causal effects can vary between individuals, which necessitates the estimation of ITEs. Treatment effects can vary between individual patients; therefore, applying a single treatment effect as CATE to all individual patients is inappropriate [47,48] as some patients may gain more or less benefit from treatments. Thus, the estimation of ITE to identify at-risk patients most likely to benefit from treatment is a major goal for stratified and precision medicine approaches. Estimating ITEs requires larger sample sizes, as individual-level estimates are less precise than aggregate-level estimates [49]. A covariate shift may result from unobserved counterfactual data but this is minimized using a weighted split-CQR approach [23].

We believe that the clinical implications of our study are significant, as understanding the causal relationships and individual treatment effects of stroke risk factors can directly influence patient care by providing more precise and personalized risk assessments. Additionally, we can conduct reviews and quality assessments of current patients in the clinic to determine who should receive further treatment. These methods enable clinicians to identify high-risk patients who would benefit most from targeted interventions, like antiplatelet therapy, thereby optimizing treatment strategies and improving patient outcomes. The use of real-world data ensures that our findings apply to everyday clinical practice.

Our study has some limitations. First, we used real-world data rather than RCT data, thus some important covariates were not previously planned, measured, and collected as part of routine clinical evaluation and were therefore unavailable for ITE estimation. Second, we acknowledge the possibility of unmeasured confounders in the observational dataset. Future studies could benefit from incorporating more granular clinical data, such as detailed medication records, laboratory results, and lifestyle factors, to mitigate potential confounding. Third, the models used for estimating ITEs were trained and validated in only a single setting, thereby limiting their generalizability. Future research should focus on validating the models in diverse settings with different patient populations or hospitals. This external validation would help to determine whether the models’ predictive performance and the estimated ITEs hold true across various contexts.

This study provides comprehensive causal estimates of AF, HT, DLP, and DM on stroke using various advanced statistical and machine learning methodologies. The consistent results across multiple analytical approaches and this study’s alignment with a standard cohort study reinforce the robustness of our findings. AF and HT emerged as significant risk factors for stroke, with DM showing a moderate effect, while DLP had minimal impact. Notably, the use of Dragonnet and conformal inference techniques allowed us to accurately estimate ITEs, highlighting that several high-risk patients who did not take antiplatelets at the time of data recorded, particularly those with DM or DM combined with HT, could potentially benefit from antiplatelet therapy. This suggests that personalized treatment strategies could be pivotal in reducing stroke risk among these patients.

The findings underscore the significance of individualized risk assessment and treatment personalization in clinical settings. Future research should focus on integrating these advanced causal inference models into routine clinical practice to enhance treatment outcomes for high-risk stroke patients. Additionally, the use of real-world data provides valuable insights but also presents challenges related to unmeasured confounding and data quality. Addressing these challenges in future studies will be crucial for advancing our understanding and improving stroke management strategies.