Comparison of B-cell depletion versus natalizumab for treatment of multiple sclerosis: A semi-supervised causal analysis

Ethics Approval

The University of Pittsburgh Institutional Review Board approved the study protocols (STUDY19080007, STUDY21030127). All registry participants provided written informed consent. EHR research protocol was deemed exempt.

Data Source

We used inpatient and outpatient EHR data from the University of Pittsburgh Medical Center (UPMC). Codified data (available from January 1, 2004 to December 31, 2022) contained demographic (e.g., age, gender, race and ethnicity) and clinical information, including daily counts of codes indicating diagnoses (e.g., International Classification of Disease [ICD]), procedures (e.g., current procedural terminology [CPT]), prescriptions (e.g., RxNorm), and laboratory tests (e.g., Logical Observation Identifiers Names and Codes [LOINC]). To organize codified data, we mapped all ICD codes to PheCodes,³³ consolidated CPT procedure codes using the Clinical Classifications Software (CCS) for Services and Procedures, and grouped electronic prescription codes to the RxNorm ingredient level (see Supplementary Materials- EHR Data Pre-processing). Narrative data (free-text clinical narrative available from January 1, 2011 to December 31, 2022) underwent a previously validated natural language processing pipeline according to unified medical language systems to generate concept unique identifiers (CUIs).^34,35

We used EHR-linked registry data of a subset of patients enrolled in a clinic-based MS cohort, Prospective Investigation of Multiple Sclerosis in the Three Rivers Region (PROMOTE, Pittsburgh, PA), between 2017 and 2023. Registry data included longitudinal outcomes (i.e., PDDS) as well as past and present DMT records, which were integrated with prescriptions from EHR to assign treatment arms.

Cohort Derivation

The cohort follow-up window (2014-2022) strikes a balance between having sufficient sample size and mitigating temporal shift in clinical practice. We used the Knowledge-Driven Online Multimodal Automated Phenotyping (KOMAP) algorithm³⁶ to identify patients with an MS diagnosis in the EHR. Among 17,827 potentially eligible patients, 3,039 pwMS received their first target treatment (BCD or NTZ) between January 1, 2014, and December 31, 2022. We excluded 993 patients who had received other highly effective DMTs (e.g., alemtuzumab, cladribine, mitoxantrone) or who received target treatment prior to January 1, 2014, as well as 76 patients who ever received chemotherapies (e.g., cyclophosphamide, methotrexate). We further excluded 217 patients who switched to another DMT after target treatment initiation.

Finally, we excluded 15 patients with missing demographic or clinical features. This final cohort included 1,738 MS patients, including 667 registry-enrolled participants.

Treatment Assignment

We compared two highly effective DMT mechanistic classes: BCD (i.e., ocrelizumab, ofatumumab, rituximab) and NTZ. While rituximab is not officially approved for MS, it has been widely used in clinical practice. Ublituximab was approved after data extraction. We assigned patients in the BCD or NTZ treatment arm, if they first received either target treatment between January 1, 2014 and December 31, 2022. The date of first receiving BCD or NTZ was the target treatment initiation date.

Disability Outcomes

To assess treatment effectiveness, we used sustained changes in disability based on PDDS^37–39. As a validated scale of patient-reported disability in MS, PDDS scores range from 0 to 8, with higher scores indicating greater disability (of gait impairment). Due to sparse observations of observed PDDS scores of 8 in the dataset (<1%), we combined PDDS scores of 7 and 8 into a single category for analysis.

First, we categorized patients into one of four groups: “sustained worsening” (≥1 PDDS score increase following treatment initiation, sustained for at least 6 months), “sustained improvement” (similarly defined but for ≥1 PDDS score decrease), “no sustained change”, and “insufficient information”. Specifically, adjudication of sustained changes required at least three PDDS scores, separated by at least 6 months between two consecutive scores, within 3 years after treatment initiation. Sustained worsening or improvement required two sequential PDDS scores that were either greater or lesser, respectively, compared to the first observed PDDS score after treatment initiation. When a patient experienced both sustained worsening and sustained improvement within the study period, we considered the first occurrence of sustained change as the outcome (i.e., worsening or improvement, whichever occurred earlier). We refer to patients with “insufficient information” as “unlabeled outcome” or as having “no labeled outcome”, and the other three categories as having “labeled” outcomes. Labeled outcome data are available only in a subset of registry-enrolled participants.

Next, we assessed sustained worsening and sustained improvement separately. When evaluating sustained worsening, we merged the other two labeled categories of patients who did not experience sustained worsening as “no sustained worsening”: i.e., patients with “sustained improvement” and “no sustained change”. When evaluating sustained improvement, we likewise merged the other two labeled categories as “no sustained improvement”: i.e., patients with “sustained worsening” and “no sustained change”.

Confounders

To control potential confounding, we adjusted for baseline features based on standard clinical trial practice, healthcare utilization, and knowledge-derived features from EHR. We obtained the covariates primarily from EHR and additionally from registry data, which provided a complementary source for registry-enrolled patients.

Standard clinical trial features included demographics (i.e., age at target treatment initiation, self-reported gender, race and ethnicity), disease duration (i.e., from the first MS PheCode to target treatment initiation), follow-up duration (i.e., from the first clinical encounter to treatment initiation), prior DMT use duration (i.e., from the earliest record of prior DMT use to treatment initiation; setting as 0 for those without any prior DMT).

As baseline disability status at treatment initiation was infrequently available, we calculated a baseline PDDS risk covariate for all patients using a PDDS imputation model and each patient’s available clinical history prior to target treatment initiation. PDDS imputation was developed using a separate dataset of PDDS scores unused in this causal analysis. See Supplementary Materials-PDDS Imputation Model, Supplementary Tables 1a-b, 2; Supplementary Figures 1a-b) for details on the PDDS imputation method and prediction performance.

Healthcare utilization features indicate not only overall health status but also potentially missing information in the fragmented US healthcare landscape. To capture total baseline healthcare utilization, we included a covariate of the counts of all observed EHR (codified and narrative) features during distinct clinical encounters that occurred prior to target treatment initiation. Further, we included a covariate of the counts of distinct clinical encounters with an MS PheCode (based on distinct code-date pairings, i.e., the number of unique dates with MS PheCode).

To select EHR features informative of MS and disability status while reducing feature dimensionality, we used an online narrative and codified feature search engine (ONCE), powered by pre-trained multi-source knowledge graph representation learning of EHR concepts.³⁶ See detailed methods in Supplementary Materials-Knowledge Graph-Derived Feature Selection and the list of included EHR features in Supplementary Tables 3a-d. To address feature sparsity, we excluded rare features with <10% frequency in the cohort. We used normalized aggregated counts of diagnoses (PheCode), procedures (CCS codes), prescriptions (RxNorm), laboratory tests (LOINC), and clinical narrative features (CUIs) before treatment initiation. Normalization was achieved by dividing the total occurrence of each feature by the total health care utilization.

Statistical Analysis

To conduct causal analyses, we employed a doubly-robust, semi-supervised estimation method.^40,41 To emulate RCT using observational data, it was crucial to adjust for confounders associated with both treatment assignment (i.e., BCD vs NTZ) and outcome (i.e., sustained change in PDDS). With confounder adjustment, the doubly-robust method yields valid estimates despite misspecification in either the propensity score model (i.e., treatment assignment given baseline confounders) or the outcome model (i.e., outcome given baseline confounders). In contrast to the standard doubly-robust methods that rely solely on labeled data (i.e., observed outcome), this semi-supervised method incorporates both labeled and no labeled outcome observations to improve efficiency in estimating comparative treatment effects. An overview of the study design is presented in Figure 1.

• Download figure
• Open in new tab

Figure 1.

Study design. (A) Schematic overview. (B) Flowchart of cohort derivation steps (boxes with blue outline) and the modeling steps in the causal analysis (green boxes). The “Non-Registry - MS” panel represents patients with a KOMAP algorithm-phenotyped MS diagnosis who were not enrolled in the clinic- based MS registry. The “Registry - MS” panel represents patients with neurologist-confirmed MS diagnosis who were enrolled in the clinic-based MS registry (with linked EHR data). Sample sizes for the “Received” DMT categories correspond to patients with a first record of BCD or NTZ between January 1, 2014 and December 31, 2022. See Methods for details on the cohort derivation. “Labeled outcome observations” correspond to patients with “sustained worsening”, “sustained improvement”, or “no sustained change” in PDDS score, while “no labeled outcome observations” correspond to patients with “insufficient information” for assessing sustained changes. Abbreviations: BCD, B-Cell Depletion therapy. DMT, Disease Modifying Therapy. EHR, Electronic Health Record. NTZ, natalizumab. PDDS, Patient Determined Disease Steps.

The semi-supervised estimation method consists of three steps. First, we fit initial propensity and outcome regression adaptive Least Absolute Shrinkage and Selection Operator (LASSO) models to calculate a modified version of a treatment propensity weight, specifically a kernel- smoothed propensity score, i.e. the double-index propensity score (DiPS). Second, we developed an outcome imputation model that predicts sustained-change (in PDDS) outcomes for data with no labeled outcome, which were then incorporated into the final doubly-robust estimation. This outcome imputation model can use additional, surrogate information available after target treatment initiation. Finally, we estimate the treatment effect using a doubly-robust inverse propensity weighting (IPW) estimator, using the DiPS as balancing weights and incorporating imputed outcomes from the second step.^40,41 This doubly-robust IPW estimator ensures consistency under correct specification of either an initial propensity score model or initial outcome model. Additionally, it is robust to misspecification of the outcome imputation model and generates consistent estimates as if the true outcomes were observed.

In the second step above, the outcome imputation model of sustained change used updated covariates including post-target treatment initiation information. This imputation model included the same set of expert-defined, codified, and narrative features included in the initial outcome model, with codified and narrative features now aggregated through 1-year after target treatment initiation. We then included three additional covariates post-target treatment features, two related to target treatment group and a third year-1 PDDS risk covariate. We included a binary covariate of the target treatment group and a weight of the inverse DiPS of a patient’s assigned target treatment group. For patients with observed PDDS scores within 1 year post- treatment initiation, the year-1 PDDS risk was the mean of observed PDDS scores. For patients without PDDS scores within the same window, the year-1 PDDS risk was the predicted PDDS score using an independently developed PDDS imputation model using codified and narrative feature data through 1-year post-target treatment. We then trained a LASSO logistic regression model of sustained change using labeled observations. We then used this model to impute the sustained change outcome when there was no labeled outcome.

The final comparative treatment effect estimation between BCD and NTZ was based on DiPS- weighted outcomes on all patients, used the observed outcome for labeled data and imputed outcomes for observations with no labeled outcome. We employed a perturbation resampling procedure to construct confidence intervals for estimates and calculated two-sided p-values through confidence interval inversion.^41,42 We considered p-values below 0.05 as statistically significant after applying a Holm correction. Further, we conducted a power evaluation based on a normal approximation, using the standard error estimated from perturbation analysis. We compared the results and efficiency of the semi-supervised analysis with an observed-case analysis (i.e., using only labeled data) using the final, doubly-robust IPW estimator. See details in Supplementary Materials-Detailed Causal Analysis.

Sensitivity Analyses

The pre-processed data for the main analyses used all available pre-treatment data to calculate pre-treatment features, trimmed rare EHR features with less than 10% frequency, chose the 6-months window between post-treatment PDDS scores to ascertain sustainment of disability change, and used knowledge graph for initial feature screening. To validate the main findings, we performed a series of sensitivity analyses by altering each of these parameters. Specifically, we repeated the analyses using 6-month and 12-month pre-treatment periods, applying a less stringent 5% frequency threshold of EHR features, categorizing sustained-change outcome status with a 3-month window between PDDS scores, and using all EHR features rather than only knowledge graph-derived features. Finally, we repeated the main and sensitivity analyses to sustained changes within 2 years (rather than 3 years) after treatment initiation.

Data Availability

The code for analysis and figures is available at <https://github.com/domdisanto/BCD_NTZ_SemiSupervisedCausal>. De-identified data are available upon request to the corresponding author and with permission from the participating institutions.

Patient Characteristics

When assessing the cohort characteristics (Table 1), registry and non-registry patients are comparable except for healthcare utilization, where registry-enrolled patients demonstrate higher mean total baseline healthcare utilization. In the overall study cohort, 1,245 (71.63%) patients received BCD and 1,478 (85.04%) did not have labeled sustained-change outcomes. Figure 2 displays the relative average or proportion of characteristics by treatment class and by sustained change outcome category as compared to the overall mean. When compared to patients in the NTZ group, patients in the BCD group had longer disease duration (5.62 vs 4.33 years), follow-up duration (9.82 vs 8.05 years), higher mean total baseline healthcare utilization (1,561 vs 900 distinct codes and CUIs), higher baseline PDDS risk (2.13 vs 1.88), and fewer women (68.35% vs 80.93%). Patients with sustained disability improvement had lower baseline and year-1 PDDS risk and were more often women. Additional summaries, stratified by outcome category and target treatment group, are available in Supplementary Tables 4a-c. Baseline covariates are similarly compared after balancing by DiPS weights in Supplementary Tables 5a-b, demonstrating appropriate covariate balance by target treatment class.

• Download figure
• Open in new tab

Figure 2.

Cohort characteristics by treatment class and outcome category. Each square reports the difference between the column-wise subgroup mean and the overall cohort mean for the variable listed in each row. Red squares indicate values for a sub-group being greater than the overall mean, while blue squares indicate values being less than the overall mean. The left panel summarizes differences across the four sustained change categories. The right panel summarizes differences across treatment assignments. For categorical variables (Women, White, Non-Hispanic race/ethnicity), the mean corresponds to the percentage of the category. Abbreviations: mo, month. PDDS, Patient Determined Disease Steps. BCD, B-Cell Depleting drugs. NTZ, Natalizumab.

Comparative Effectiveness Analysis

The main comparative effectiveness analysis results are presented in Table 2. Higher PDDS scores indicate greater patient-reported disability. With NTZ as the reference treatment class, a positive ATE for the sustained worsening outcome would indicate BCD as less effective than NTZ in reducing patient-reported disability worsening. Conversely, a positive ATE for the sustained improvement outcome would indicate BCD as more effective than NTZ in promoting patient-reported disability improvement. We observed no statistically significant difference between BCD and NTZ for either sustained worsening or sustained improvement during the 3- year post-treatment evaluation.

Features associated with treatment assignment and sustained-change in disability outcome informing the final, doubly-robust IPW estimator in the main analysis are shown in Figure 3. The same propensity score model and thus the informative features in the DiPS were used in analysis of both sustained-change outcomes. In the treatment model, men (as self-reported gender), higher baseline total healthcare utilization, mental disorders (PheCode: 306) and older age at target treatment initiation were most associated with a higher probability of receiving BCD than NTZ. Features associated with increased risk of sustained disability worsening included muscle spasms (PheCode 350.1), dysuria (PheCode 599.3), distress (C0231303), and leukoencephalopathy (C0270612). Few features were selected in models of sustained disability improvement with the largest coefficients observed for treatment with NTZ and metabolic diagnosis procedures (C4263342).

• Download figure
• Open in new tab

Figure 3.

Coefficients in the initial propensity score model (A) and two initial outcome models (B) fitted to inform the kernel-smoothed propensity score estimator. The same treatment model was used for final causal analysis of each outcome. These models collectively informed the final inverse probability weight that was used to weight/balance for the final ATE calculation, but no balancing/weighting was used to calculate the coefficients shown here. Only variables with non- zero coefficients are included in each panel. Dot size represents the magnitude of the coefficient, while color corresponds to both effect size and direction. In Panel A, covariates with larger (positive) coefficients are associated with greater probability of receiving BCD than NTZ while the inverse is for smaller coefficients. Detailed descriptions of CUI, PheCode, and RxNorm codes are included in Supplementary Tables 3a-d. All covariates were standardized (to mean 0, unit variance) prior to model fitting. Abbreviations: CUI, Concept Unique Identifier.

Sensitivity Analyses and Power Evaluation

Sensitivity analyses reducing the pre-treatment periods, lowering the frequency threshold of EHR features, reducing the interval between PDDS scores in the sustained-change outcome definition, and using all available (rather than only knowledge graph-derived) EHR features produced treatment effect estimates and confidence intervals similar to the primary results (Supplementary Figures 3a-b, Supplementary Tables 6a-b), as did analyses using the 2- year post-treatment evaluation window (Supplementary Table 7, Supplementary Figure 2).

To quantify the advantage of the semi-supervised method leveraging observations with no labeled outcome, we compared the main semi-supervised analysis to an observed-case analysis (i.e., using labeled data only) in the final ATE estimator (Supplementary Tables 8a-b). Fixing power at 80%, we reported the minimum detectable effect size for the risk difference in sustained improvement and sustained worsening between BCD and NTZ. Notably, the semi- supervised estimate is powered to detect a risk difference half the effect size of the observed- case estimate, indicating a meaningful gain in power when using the semi-supervised estimator.