In this study, we utilized the claims dataset obtained from the National Health Insurance Service (NHIS) in South Korea, which ensures comprehensive healthcare coverage and collects health records for over 99.6% of the Korean population (Ro et al., 2019). These anonymized records, available for research, include demographic details, diagnoses, prescriptions, and health screening data of individuals.

For discovery and sensitivity analyses, we examined two specific cohorts: the customized cohort (CC) and the national sample cohort (NSC), both sourced from the NHIS. The CC, derived from the pre-established NHIS data based on researcher-specified criteria (Kyoung and Kim, 2022), provides unrestricted access to clinical and demographic data but limits prescription information to only the drugs requested. Due to its larger and more representative sample size, the CC was used for discovery analysis. The NSC, a randomly selected subset representing 2% of the national population (Lee et al., 2017; Kyoung and Kim, 2022), contains comprehensive records for all medications (Lee et al., 2017). We used the NSC for the sensitivity analysis to assess whether the results from the discovery analysis remained consistent when additional adjustments for concomitant medications were made.

We focused on patients diagnosed with OA, classified under International Classification of Diseases (ICD) codes M15–M19 from 2010 to 2014. The criteria for inclusion were as follows: no prior OA diagnosis from 2002 to 2009, age over 50 years, availability of health screening data without missing values, and no history of rheumatoid arthritis (ICD codes M05, M06), infectious or inflammatory arthritis (ICD codes M00–03, M10-14), traumatic arthritis or prior fractures (ICD codes S82, M17.2, M17.3), or osteonecrosis (ICD code M87).

We categorized AEs using the first three characters of ICD codes, with each unique three-character ICD code treated as a distinct potential AE related to NSAID use under the following conditions.

1. We excluded codes unrelated to our main exposure, as listed in Supplementary Table S2.

3-1. For each ICD code, we calculated the prevalence of NSAID prescription on the same date conditional on the diagnoses with the ICD, and it was compared with the marginal prevalence of NSAID prescription. If the ICD code is an indication for NSAID use, the former is expected to be larger than the latter. The conditional probability of drug prescription given a specific diagnosis code is calculated as follows: P NSAID prescription | ICD diagnosis = N umber   of   NSAID   prescriptions   issued   on   the   day   t h e   ICD   code   was   recorded N umber   of   t i m e s   t h e   ICD   codes   were   recorded

The latter during the study period was 16.16%, and it was calculated as P NSAID prescription = N umber   of   NSAID   prescriptions   occurring   on   the   same   day N umber   of   overall   hospital   visits

Chi-square test was used to compare it with P NSAID prescription | ICD diagnosis at the 0.05 significance level. Consequently, prescription rates for 83 ICD codes exceeded the average (Supplementary Table S3).

3-2. If certain ICD codes are indications for NSAID use, it is expected that the number of prescriptions increased significantly after a recorded diagnosis. For each ICD, we compared its prevalence 10 days before an NSAID prescription with its prevalence 10 days after NSAID prescription by Poisson test. This comparison led to the identification of 183 codes (Supplementary Table S4), which were subsequently removed from the analyses. After these filters, 246 three-digit ICD codes remained and were used as potential AEs for our statistical analyses.

For each potential AE, its occurrence was detected using the corresponding ICD codes. We considered events occurring 1 year after the index date to prevent immortal time bias. In Korea, NSAIDs are usually co-prescribed with gastroprotective medicine to prevent gastrointestinal (GI) complications (Jeong et al., 2023). These prescriptions should be accompanied by GI-related ICD codes such as K20–K31, and GI events occurring on the same day as an NSAID prescription were not considered AEs. After these filters, 246 three-digit ICD codes remained and were used as potential AEs for our statistical analyses.

We used multiple analyses to screen for AEs. For discovery purposes, we used CC data from 888,909 newly diagnosed (2010-2014) OA patients, employing both cohort and case-only study designs. The main exposure was NSAID prescriptions associated with OA, considering various subtypes classified according to the anatomical therapeutic chemical (ATC) classification system (detailed in Supplementary Table S5). The observation period ended in 2018.

In the cohort design, each participant prescribed NSAIDs was matched with five non-NSAID users based on propensity score matching with age, sex, and body mass index (BMI). The index date was set at the issuance of the first NSAID prescription following an OA diagnosis; for matched non-users, this date aligned with their medicated counterparts. Collectively, we included 120,691 NSAID users and 548,978 non-users. Exposed groups were defined using medication possession ratios (MPRs), which are used to quantify the proportion of days a medication was prescribed annually, serving as a measure of drug compliance. Patients were categorized into non-medication groups (MPR <0.1 throughout the observation period) and medication groups (MPR ≥0.1 at any point during the study).

In the case-only design, all participants (N = 888,909) were included without segregation into medication/control groups, with the observation start 1 year before the OA diagnosis to avoid potential early clustering of AEs.

For sensitivity analyses, significant results from the discovery dataset (CC) were replicated after adjusting concurrent treatments using a separate NSC dataset. We expanded the study to include patients newly diagnosed with OA from 2010 to 2018 to ensure a large sample size, using the same criteria and index dates as the original dataset, ultimately involving 23,138 medicated and 60,669 non-medicated patients.

To adjust for confounding effects in cohort design, several potential confounders such as age, sex, BMI, Charlson comorbidity index (CCI), income, and residence were incorporated as covariates. Specifically, age, BMI, and CCI were included as continuous-scale covariates, whereas sex, income, and residence were considered as categorical variables. Income was divided into 20 groups, and residence was classified as a metropolitan city, small city, or rural area. These variables were adjusted for at the ordinal and nominal scales, respectively.

Case-only design inherently accounts for covariates within individuals and, thus, does not require additional adjustments. The OA diagnosis date was used to account for changes in the time-varying condition related to OA itself.

In the sensitivity analyses, we used the same design and covariates as the cohort study, with the additional adjustment for concomitant drugs related to the classified results.

For demographic statistics, continuous variables were tested using t-tests, and categorical variables were tested using Chi-square tests.

For the cohort study, we conducted Cox proportional hazard analyses with the onset age of potential AEs as the outcome. For each potential AE, a 1-year washout period was applied, and participants with any event during that period were further filtered out on a per-AE basis. Annual MPRs were treated as time-dependent exposures, leading to the execution of an extended Cox proportional hazards model. The effect of NSAID on each AE was estimated after adjusting for covariates. Subsequently, we identified an association between each identified AE and the five most commonly prescribed NSAIDs: aceclofenac, meloxicam, loxoprofen, celecoxib, and naproxen. We used the same Cox proportional hazard model, with the primary exposure variable shifted to these five types of NSAIDs. Subsequently, pairwise comparison was performed using a one-sided Wald test.

For the case-only study, self-controlled case series (SCCS) analyses were performed to explore the significant association of each potential AE (Ahn et al., 2023). Similar to the cohort study, a 1-year washout period was implemented. In the SCCS framework, the risk period was defined from the initiation of NSAID use, extending to 6 months after cessation of the medication. Within this timeframe, “current use” was defined as the actual duration of medication use, with an additional 14-day period designated as “recent use.” The remaining days within the risk period were classified as “past use.” To address potential biases arising from reduced observed periods, particularly those resulting from high fatality rates of severe AEs, an extended SCCS model was employed in relevant cases. The ICD codes consistently identified in both analyses were defined as AEs of NSAIDs.

A Bonferroni-adjusted significance level of 0.05 was established. Given the large sample size, we also set an additional threshold of |beta| > 0.2 to ensure substantial differences. These criteria were consistently applied to both cohort and case-only studies, and only results meeting them were considered significant.

In the sensitivity analyses, we first conducted a replication using the same Cox proportional model and covariates as in the cohort study from the discovery analysis. We Subsequently, we added drug-related covariates to check if the results remained consistent. Due to the insufficient sample size to perform a case-only design, sensitivity analyses were focused on the cohort study. For clarity, all references to “cohort study” in the following text will specifically refer to the cohort study in the discovery analysis, while “sensitivity analysis” will be used as described. A detailed study flow and diagram are presented in Figure 1 and Supplementary Figure S1.

The statistical significance level was set to 0.05. All analyses were performed using R software (version 4.0; R Project, Vienna, Austria), SAS Enterprise Guide (version 7.15; SAS Institute Inc., Cary, NC, United States), and Rex (version 3.5.3; RexSoft, Seoul, Republic of Korea).

The study population characteristics are summarized in Table 1. The mean age of the medicated groups was significantly higher than that of the non-medicated groups (P < 0.001). Despite significant differences across all variables, standardized mean differences for matching variables (age, sex, and BMI) were confirmed to be <0.2, ensuring a balance between the groups. Additionally, the frequencies of subtypes and daily doses of prescribed NSAIDs are detailed in Supplementary Figure S2. Aceclofenac was the most prescribed NSAID for OA, followed by meloxicam, loxoprofen, celecoxib, and naproxen.

We considered 246 potential outcomes in discovery analyses, with cohort study results for all included ICD codes presented in Supplementary Figure S3. Eight codes demonstrated significance at the 0.0002 level (equivalent to the 0.05 BF-adjusted significance level), consistently showing the same direction of effect. Detailed information on significant AEs is provided in Figure 2. These AEs were categorized based on similarities in diagnosis: anemia (D50, D62, and D64), cardiac effects (I10 and I50), cerebrovascular disease (I63), and renal toxicity (E87 and N17).

The highest adjusted HR (aHR) derived from the cohort study for the general NSAID MPR was observed for acute posthemorrhagic anemia (D62; adjusted HR [aHR] = 2.45, 95% confidence interval [CI] = 1.71–3.52) and acute renal failure (ARF) (N17; aHR = 2.16, 95% CI = 1.68–2.77) (Figure 2A). Consistent with the cohort study, acute posthemorrhagic anemia and ARF exhibited the highest relative incidences (RIs) in the case-only study. The RIs for these diagnoses were significantly higher than one during both current and recent periods, indicating that the incidence of AEs lasted up to 2 weeks post-medication. Elevated RIs became not significant for the past period, from 2 weeks to 6 months post-NSAID intake, for all AEs. (Figure 2B). Additionally, the sensitivity analysis adjusting for concomitant medication did not alter the significance of the results (Figure 2C). The list of concomitant drugs is presented in Supplementary Table S6.

Figure 3 displays the estimated effects of individual NSAIDs, and Supplementary Figure S4 presents the results of the pairwise comparisons. Figure 3A shows the outcomes from the cohort analyses. For anemias, non-selective NSAIDs exhibited significantly higher risks than celecoxib (pairwise P < 0.05). When examining hypertension (I10), the risks associated with all NSAID subtypes were comparable (pairwise P > 0.05), with loxoprofen demonstrating a relatively lower hazard for heart failure (HF) (I50; pairwise P < 0.05). Celecoxib was the only NSAID that showed a significant association with cerebrovascular disease, although the risk was not significantly higher than that of other NSAIDs (pairwise P > 0.05). In terms of renal toxicity, most NSAIDs demonstrated significance, except for meloxicam. Notably, naproxen exhibited a substantially high HR (aHR 4.7, 95% CI 2.16–10.24) (P < 0.0001) (Figure 3A; Supplementary Figure S5); however, the risk was comparable with those of celecoxib and loxoprofen (pairwise P > 0.05). The results from the case-only studies were consistent with those from the cohort ones (Figure 3B).