Incident LC patients and matched controls were selected from the National Health Insurance Service-National Sample Cohort (NHIS-NSC) of South Korea. The NHIS established an administrative database for research purposes, i.e., the National Health Information Database (NHID), which stores all healthcare service records from the general population [15,16]. From the NHID, the NHIS constructs and provides researchers the NHIS-NSC, a representative administrative dataset, for guiding health policy and carrying out biomedical research. A detailed cohort profile was published by the NHIS [16]. The NHIS-NSC consists of one million people selected by a systematic stratified random sampling method from a target population of 48,222,537 in 2006. Under the aegis of a compulsory social insurance system instituted by the National Health Insurance Act, all subjects could be followed until 2015, except in cases of disqualification due to death or emigration. Data from 2002 to 2005 were included to provide additional past medical information. As a single payer, the NHIS maintains all personal information, demographics, and medical treatment data for the entire Korean population. Key variables in the NHIS-Senior include all inpatient and outpatient medical claims data, such as treatment procedure codes, prescription codes, and diagnostic codes. This study protocol was approved by the Institutional Review Board of Chungnam National University Hospital (IRB No.: 2019-10-001).

The incident LC cohort was constructed from the NHIS-NSC base cohort. To be enrolled, subjects must fit the operational definition of new-onset LC, i.e., first-time admission to an acute care hospital or two or more outpatient visits with primary or the first secondary diagnostic codes of LC (International Statistical Classification of Diseases and Related Health Problems, 10th Revision [ICD-10] K702, K703, K717, K74, K761, K765, K766). The decompensated LC cohort comprised those with K702, K703, K717, K74, K761, K765, K766 codes accompanying I85.0, I86.4, I98.3, K72, K72.1, K72.9, K76.7, R17, R18 codes, and the alcoholic cirrhosis cohort was defined as patients with simultaneous codes of K70.2 and K70.3. For definition of the ICD-10 code for LC patient groups, validation studies and a South Korean study were referenced [17-20]. Among these patients, some were excluded to maintain the study design validity. To select new-onset LC cases, patients having medical claims of LC between 2002 and 2006 (a 5-year washout period) were considered to have a diagnosis of LC previously, and were excluded. The date of LC incidence was defined as that of the first admission or the second outpatient visit, whichever occurred first. The second outpatient visit was considered the date of diagnosis of LC to avoid immortal time bias [21], which is the interval between the first and the second outpatient visit. LC patients younger than 18 years at the time of LC incidence were also excluded. From January 1, 2007, new-onset LC patients were enrolled in the incident LC cohort at the index date (time zero) and followed up.

Information on death from Statistics Korea was individually linked using unique personal identification numbers [16,22]. Under a de-identification policy, the calendar year and month of death were provided by the NHIS, after masking the exact date. By law, all deaths are required to be reported to Statistics Korea. Suicide was identified by the set of deaths assigned ‘Intentional self-harm’ (ICD-10 codes, X60 to X84).

Although the NHIS-NSC was constructed retrospectively, the design of this study was planned to mimic a prospective study. First, a time-dependent propensity score was calculated, followed by risk-set matching [23,24]. To adjust for confounding effects, the association between LC and suicide was examined by time-dependent propensity score matching [24]. Hazard components (as propensity score) were estimated using Cox’s proportional hazards model with January 1, 2007 as baseline and LC as event. All variables included in Table 1 were independent variables and were collected from 3 years (2004 to 2006) before baseline. Age and square of age were included as continuous variables, whereas sex, household income level (decile), registered disability, residential district as fixed effect, Charlson Comorbidity Index (CCI), numbers of hospital admissions and outpatient visits, and past medical history were included as categorical variables. Each subject’s number of comorbidities was assessed by diagnostic codes using the Quan ICD-10 coding algorithm of CCI score [25]. The presence of CCI disease-constituting categories and the presence of a past medical history of malignant neoplasm, congestive heart failure, chronic obstructive pulmonary disease, or alcoholic liver disease were defined by at least two outpatient visits or one admission upon primary or first secondary diagnosis. Prescription over 90 days for antihypertensive, antidiabetic, or lipid-lowering agents was accounted for as patients using corresponding medications.

Second, to emulate a prospective study, the patients were matched to subjects at risk of LC by age, sex, and use of anti-depressants at the calendar date (time zero) when the first LC patient was identified. This method of risk-set matching was repeated for subsequent patients until the last LC patient [23,26,27].

Third, 1:2 matching on propensity score was performed sequentially for each risk set using a nearest neighbor matching algorithm with a maximum caliber of 0.1 of the hazard components. To allow matching to be independent of future events, the matched-control subjects could be either those who never developed LC or those who had not yet developed LC. Thus, an LC patient in the incident LC cohort could enter the study as an LC patient or a matched control for another LC patient whose time zero was prior to that of the initial patient [28]. Next, to yield non-overlapping samples from the risk set, the matched subjects were removed from the next risk sets. The process was performed iteratively with the next risk sets until there were no more treated patients in the risk set.

Using the final matched cohorts, statistical tests for the association between incident LC and risk of suicide were performed considering the statistical nature of matched-pair analysis. To assess covariate balance between treatment groups, baseline characteristics were compared with standardized differences where a difference less than 0.1 (10%) is generally considered negligible [29,30]. The cumulative incidence curve of suicide was obtained by Kaplan-Meier survival curves and a stratified log-rank test was used to compare the Kaplan-Meier curves of the matched cohort [31]. Cumulative incidence of suicide and 95% confidence interval (CI) was calculated by the product limit (Kaplan-Meier) method of survival probability. Incidence rate (IR) of suicide and 95% CI were calculated based on a generalized estimating equation with a Poisson distribution and expressed as number of suicides per 100,000 person-years. The effect size was presented as hazard ratio (HR) using Cox’s proportional hazards model with a robust variance estimator that accounts for clustering within matched pairs [30,31]. Time zero was set to the date of LC incidence for both LC patients and their matched controls. Survival time used in analyses was defined as months from time zero to time of death or December 31, 2015, whichever came first. Proportional hazards assumption was assessed by graphically using log transformation of the negative log of estimated survivor function, time-dependent explanatory variables, Schoenfeld residuals, cumulative sums of martingale residuals, and a supremum test. In addition to these cause-specific models for competing risks, a Fine and Gray subdistribution hazard model was performed as a sensitivity analysis. Sub-group analyses within categorized age groups (18–49, 50–64, and 65 years or more), sex, and follow-up time frames (0–24, 25–60, 61 months or over) were performed. To compare the suicide rate according to severity of cirrhosis, the analysis was performed by dividing patients into matched cohort, compensated LC, and decompensated LC. Moreover, to additionally perform stratified analysis on patients with depression, as a coexisting mental illness, the suicide rate according to presence or absence of LC was analyzed by categorizing whether anti-depressant was administered. The statistical analyses were conducted using SAS Enterprise Guide version 7.1 (SAS Institute, Cary, NC, USA). P-values <0.05 were considered statistically significant.

From January 1, 2002 to December 31, 2015, a total of 10,594 patients met the inclusion criteria. Among these patients, LC patients during the washout period of 2002 to 2006 (4,354 patients) were excluded. During risk-set matching, five LC patients were not matched to control subjects and were excluded, and another 331 patients entered into the study as control subjects of another LC patient were excluded.

Additionally, 87 LC patients were excluded due to a less than 1:2 matching ratio, and eight LC patients younger than 18 years were excluded. Finally, 5,809 incident LC patients and their 11,618 matched controls were analyzed in this study (Fig. 1). The mean period of follow-up was 4.20 years, generating 73,261 person years. During follow-up, a total of 61 suicide cases was identified. Table 1 shows the baseline characteristics of the matched cohorts. Mean age was 56.3 years, and 28.6% of the patients were female. The highest standardized difference in number of outpatient visits was 5.19%.

Cumulative incidence of suicide during the entire follow-up period showed a statistically significant difference between LC and matched-control cohort (P=0.002 for stratified log-rank test) (Fig. 2). The 5-year cumulative incidence risk of suicide was 0.78% (95% CI, 0.52–1.13) for the LC cohort and 0.33% (95% CI, 0.22–0.47) for the matched-control cohort (Table 2). During the entire follow-up period, among 5,809 LC patients 30 suicides were identified during 20,930 person-years (IR, 143.3 per 100,000 person-years; 95% CI, 100.2–205.1). The LC patients were 2.37 times more likely to commit suicide than their matched controls (HR, 2.37; 95% CI, 1.44–3.88). Moreover, increased suicide risk was predominantly observed within the first 2 years of follow-up (HR, 2.59; 95% CI, 1.20–5.60), among patients in the 18–49-year-old group (HR, 3.72; 95% CI, 1.45–9.56) (Fig. 3).

In subgroup analyses, non-alcoholic LC patients had a 2.04 times higher suicide risk than the matched cohort (HR, 2.04; 95% CI, 1.16–3.58), and alcoholic LC patients had a 3.29 times higher suicide risk compared with the matched cohort (HR, 3.29; 95% CI, 1.66–6.52). However, when comparing the suicide risk of nonalcoholic LC patients and alcoholic LC patients, the suicide risk of alcoholic LC patients was 1.62 times higher, although there was no statistical significance (HR, 1.62; 95% CI, 0.77–3.40; P=0.206). Moreover, compensated LC (HR, 2.03; 95% CI, 1.18–3.49) and decompensated (HR, 4.42; 95% CI, 2.05–9.52) LC patients were more likely to commit suicide compared to their matched controls. Interestingly, as a result of a comparative analysis of suicide risk between compensated and decompensated LC patients, the suicide risk of decompensated LC patients was 2.18 times higher, although the difference was not significant (HR, 2.18; 95% CI, 0.96–4.94; P=0.062). Among ever-users of antidepressants, increased risk of suicide was observed in LC patients, although the association was not significant (HR, 1.97; 95% CI, 0.55–7.13; P=0.301). Among patients not taking anti-depressants, the risk of suicide in LC patients was significantly 2.44 times higher compared to that of the matched cohort (HR, 2.44; 95% CI, 1.42–4.17; P=0.001) (Fig. 4). Results of Fine and Gray sub-distribution hazards model were similar to those of the cause-specific model.