2.1 Simulation Generation Process

We conducted a series of simulation studies to investigate the biasing effects of various chronic medication-use patterns in case-crossover and CTC designs. For each medication-use pattern, we simulated 1000 repetitions of N = 10,000 patient cohorts, where each patient had 1460 daily observations. We then restricted the analytic period to days 366–1095 to avoid unstable prevalence of exposure during the observation period [14] and designated the first year of the trimmed observation timeframe as the look-back window to be consistent with the continuous enrollment criterion needed for case-crossover studies. A schematic depiction of the simulation generation process is illustrated in Electronic Supplementary Material (ESM) Fig. 1. We applied two distinct approaches to simulate chronic medication patterns: (1) separate patterns of common chronic medication use and (2) a cohort with a combination of 15 medication patterns reflective of distribution of patterns of chronic medication use in the real world.

Fig. 1

Study design diagram for the main and the sensitivity analyses

2.2 Separate Patterns of Common Chronic Medication Use

To fully understand the impact of specific patterns of medication use, we simulated seven patterns (P1–P7) that represented the shape of five common trajectories reported in the literature for chronic medications [10]. These patterns were characterized as (1) early discontinuers, (2) rapidly declining, (3) U-shaped, (4) slowly declining, and (5) consistent adherence [10]. For each common trajectory, we generated annual medication-use patterns based on 30-day increments to represent the most common dispensing of non-transient medications with 30 day supplies [15]. Thus, each pattern could be described with a 12-digit combination of months with exposure (1) and non-exposure (0) to reflect chronic medication use. Details for each pattern and how these patterns were selected are included in ESM Table 1.

2.3 Real-World Chronic Medication Patterns

To capture real-world chronic medication use patterns, we used MarketScan data to identify real-world chronic medication patterns among SGLT2 inhibitor users. The MarketScan data is an administrative claims database in the USA that contains data including outpatient drug dispensing, healthcare utilization, patient demographics, and inpatient/outpatient claims [16]. Detailed information on identification of SGLT2 inhibitor users is presented in ESM Fig. 2. Among SGLT2 initiators from October 2015 to December 2018, we determined annual medication-use patterns based on monthly indicators of medication use. Like our simulated patterns described above, we characterized the medication-use patterns observed over each of the months as exposed (1) or unexposed (0) which resulted in a 12-digit combination of ones and zeros to represent medication use over the year. For example, if an individual was dispensed an SGLT2 inhibitor for 30 days with no subsequent dispensing, the medication pattern would be 100000000000. Among SGLT2 inhibitor users, we identified about 850 unique patterns and selected the 14 most common patterns to represent the real-world chronic medication patterns used in this approach.

In the simulated cohort, each pattern identified in the real-world data was assigned proportionally based on the observed frequency of that pattern (ESM Table 2). The 14 most prevalent patterns represented 82.5% of the population. The patterns assigned to the remaining 17.5% of patients were random patterns of medication use. To better understand each selected pattern and the expected result from this approach, we plotted the exposure prevalence of the simulated cohort and all 15 patterns (including random patterns) across 1460 days (ESM Fig. 3).

2.4 Outcome Generation

The outcome in this simulation study was assigned based on the daily probability of outcome during the event period from day 730 to day 1095. Selecting this 1-year period to assign events ensured that everyone in the analysis had at least 365 days of observation time preceding their event date. The daily probability of a binary outcome was calculated based on a function of exposure for each day (i) for person (j):

$$\text(Yij=1)= \text(\beta 0+\text\left(\beta \right)Exposure\left(ij\right))$$

The events were randomly distributed across the event period based on the daily probability. We simulated baseline risk (\(\beta 0)\) of the outcome at 0.3 and randomly selected an outcome during the evaluation window if patients were assigned multiple days with an outcome [9]. We varied the effect of exposure across a range of values (odds ratio (OR) of 0.50, 0.75, 1.00, 1.50, 2.00) in the simulations.

2.5 Statistical Analysis

Each simulation was analyzed using conditional logistic regression that compares the odds of exposure on the outcome day to the odds of exposure on a referent day 30 days prior to the outcome for the case-crossover analysis (point-in-time approach) [17]. Each of the windows were considered exposed if the drug exposure overlapped with the window. We also implemented the CTC in each simulation. We used risk-set sampling to select one control from the simulated cohort for each case. Individuals in the simulation were eligible as controls at any point before their outcome day. The OR for CTC was calculated using an interaction term between case status and exposure discordance [6]. To validate the construction of our simulated cohort, we conducted a cohort analysis that calculated the relative risk (RR) comparing exposed days to unexposed days across each of the simulations with the exception of the simulated real-world cohort.

For each analysis of a single pattern, we created a box plot of the ORs across the 1000 iterations for each pattern. We also determined the median OR and the 95% confidence interval (CI) for each of the simulations. To quantify the bias introduced by different treatment patterns, we calculated the absolute difference and relative bias between the estimated effect and the true effect across all simulated ORs. Relative bias was derived by dividing the difference between the estimate OR and the effect estimation to the true OR by the effect estimation to the true OR for each iteration ((estimate-OR)/OR). In the simulation representative of the real-world cohort, we estimated the median effect estimates and absolute difference between effect estimates and true effect and reported the 95% CI across all simulated ORs. Statistical significance was assessed at the 5% level, using a 95% CI excluding 1.00 (or a p value < 0.05).

The sensitivity analysis was conducted based on previous literature that suggested extending the risk and control window may reduce the bias for non-transient medications for case-crossover design [4, 18]. We conducted a sensitivity analysis where we extended the risk and control window from 1 day to 30 days (Fig. 1). All simulations and analyses were performed using R 4.2.2 (R Core Team 2022. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/).

2.6 Real-World Case-Example

We examined the association between SGLT2 inhibitor use and genital infections (positive outcome) and retinal detachment (negative outcome) as applied examples with chronic medication and implementation of the case-crossover design. Both outcomes were selected based on prior studies of SGLT2 medications. Previous literature has demonstrated an increased risk of genital infections among patients using SGLT2 inhibitors compared to other antidiabetic medications in both clinical trials and observational studies [12, 13, 19]. Retinal detachment was selected as a negative control outcomes in prior studies that evaluated the persistent user bias associated with the case-crossover design [8]. We used the same SGLT2 users cohort from MarketScan  to derive the real-world SGLT2 medication-use patterns (ESM Fig. 2). Patients with genital infection or retinal detachment were identified using International Classification Disease (ICD)-10 codes from inpatient and outpatient encounters. The first date of the outcome following the type 2 diabetes (T2D) diagnosis date was defined as the event date. Patients were considered exposed during a risk or control period if the day supply for an SGLT2 inhibitor dispensing extended into that period.

We used the same study designs separately for both the positive and negative outcomes that were used in the simulated approaches. First, we conducted a case-crossover design using conditional logistic regression. Second, we conducted a CTC study where risk-set sampling was used to identify non-case controls from the cohort. Each case was matched with one control based on age, sex, and cohort entry month. Once matched, controls were assigned the same event date as their matched case. The CTC estimates were calculated using the same method as the case-crossover design with an interaction term between case status and exposure discordance [6]. We implemented the same sensitivity analysis that was used in the simulations (Fig. 1). Estimates of the association between the positive and negative outcomes are presented using ORs and 95% CIs. Analyses in the applied example were performed using SAS version 9.4 (Cary, NC, USA).