Cost-effectiveness of hemodialysis and peritoneal dialysis: A national cohort study with 14 years follow-up and matched for comorbidities and propensity score

Although treatment for the dialysis population is resource intensive, a cost-effectiveness analysis comparing hemodialysis (HD) and peritoneal dialysis (PD) by matched pairs is still lacking. After matching for clinical characteristics and propensity scores, we identified 4,285 pairs of incident HD and PD patients from a Taiwanese national cohort during 1998–2010. Survival and healthcare expenditure were calculated by data of 14-year follow-up and subsequently extrapolated to lifetime estimates under the assumption of constant excess hazard. We performed a cross-sectional EQ–5D survey on 179 matched pairs of prevalent HD and PD patients of varying dialysis vintages from 12 dialysis units. The product of survival probability and the mean utility value at each time point (dialysis vintage) were summed up throughout lifetime to obtain the quality-adjusted life expectancy (QALE). The results revealed the estimated life expectancy between HD and PD were nearly equal (19.11 versus 19.08 years). The QALE’s were also similar, whereas average lifetime healthcare costs were higher in HD than PD (237,795 versus 204,442 USD) and the cost-effectiveness ratios for PD and HD were 13,681 and 16,643 USD per quality-adjusted life year, respectively. In conclusion, PD is more cost-effective than HD, of which the major determinants were the costs for the dialysis modality and its associated complications.


Methods
Establishment of the national cohort of the dialysis population. This study was approved by the ethics review board of National Cheng Kung University Hospital (A-ER-101-089) before commencement, and the methods were carried out in accordance with the approved current guidelines. The Taiwan National Health Insurance (NHI) program was initiated in March 1995, and by 2001 it provided medical services for more than 99% of the 23 million residents in Taiwan 37 . Nearly every kind of medical service can be reimbursed by the NHI, including all payments for outpatient and inpatient services, medication prescriptions and intervention procedures. Each beneficiary is required to pay a certain amount of copayment toward the utilization of the medical service. The rate of copayments ranges from less than 10% for usual outpatient visits to the highest rate of 30% for prolonged hospitalization. However, patients with specific diseases on the list of catastrophic illness, including ESRD on maintenance dialysis therapy, can be waived from copayment. Therefore, each patient with any catastrophic illness must be validated by at least two specialists to avoid abuse of the NHI program. Moreover, the waiving of copayments for ESRD under maintenance HD and PD makes our estimation of total cost comprehensive and accurate. This study was conducted using the reimbursement database of the NHI. We performed a nationwide collection of incident ESRD patients (International Classification of Diseases, Ninth Revision, Clinical Modification codes [ICD-9 codes]: 585) with certifications of catastrophic illness related to dialysis therapy (HD or PD). Patients older than 18 years and receiving dialysis for more than three consecutive months from January 1, 1998, to December 31, 2010, were enrolled in the dialysis cohort (Fig. 1). The dates of the enrollment of dialysis patients were the first date of receiving dialysis for three consecutive months. HD and PD therapies were identified according to their procedure codes. The survival of each dialysis patient was verified by cross-linkage with the database of the National Mortality Registry. Because patients with malignancy (ICD-9 codes: 140-208, 230-234) are usually associated with high mortalities, poor QOL and high subsequent costs of management for the malignancy and associated complications. In addition, there is a great variation in life expectancies (LEs), QOL and healthcare expenditures among malignancies of different organ-systems and stages. Therefore, we excluded all patients with malignancy to avoid potential confounding of estimations by cancer (Fig. 1). In addition, patients who had ever changed the dialysis mode for more than three months were also excluded because of the difficulty to quantify the effects of shifting dialysis modality at different dialysis vintages on survival, costs and QOL. Since many patients receiving renal transplantation were first treated with dialysis while waiting for an appropriate donor, they were treated as censored on the date of transplantation. However, we excluded all patients with previously performed transplantation before dialysis to avoid any potential bias in this study due to prior transplantation. Major comorbidities noted before initiation of dialysis were identified with the corresponding ICD-9 codes (listed in Supplementary Table 1), as follows: If the patient had any of these comorbidities in the hospitalization discharge codes once, or at least twice in ambulatory care with 30 days apart, within one year. Such procedures for analyzing the NHI database could improve the accuracy of measurement of comorbidities and have contributed to many high quality studies [38][39][40][41] . The end of the follow-up period was December 31, 2011. Collection of information of EQ-5D by cross-sectional sampling in prevalent dialysis patients. From February 2012 to May 2013, we conducted a cross-sectional survey of prevalent dialysis patients from 12 dialysis centers to collect the information of QOL. Informed consent was obtained from all study participants. We applied the EQ-5D to assess health-related QOL, which has five domains and can be converted into utility values for estimation of QALY 42 . Each patient was assessed by a researcher assistant who was experienced in conducting the EQ-5D survey. When necessary and applicable, information of the EQ-5D was obtained from the family members or primary caregivers of the patients as proxies. The utility value of EQ-5D ranged from 0 to 1, where 0 indicates death and 1 indicates perfect health 43 . The concomitant diseases of patients that were noted while interviewing were recorded as comorbidities. As before, patients with malignancy, dialysis mode switching or aged under 18 years old when visiting were excluded from the analysis.
Matching based on individual characteristics and propensity scores. Since the choice of dialysis modality of the physician and/or patients is based on minimizing mortality and likelihood of developing complications, as well as improving QOL, or, the overall outcomes, we first matched the two groups on major comorbidities that would result in premature mortality, and this was then followed by a matching of propensity scores to minimize potential residual confounding. The whole process thus takes the preference of physician and/or patients into consideration. A 1:1 matching of incident HD and PD patients from the national dialysis cohort was performed for age (± 2 years), sex, index year of initiation of dialysis, urbanization status, and major comorbidities, including diabetes mellitus, acute myocardial infarction, congestive heart failure, stroke, chronic liver disease, and the propensity scores up to ± 0.05. The propensity score for PD prescription was estimated with a logistic regression model 44,45 , which included hypertension, coronary artery disease, cardiac dysrhythmias, peripheral vascular disease, hyperlipidemia, and rheumatological disease as independent variables. These baseline characteristics were chosen for matching because of their association with LEs and healthcare expenditures, which were the parameters of cost-effectiveness 21,23,24,46 (Fig. 1). To control the confounding of QOL, we conducted 1:1 matching on the cross-sectional samples interviewed for QOL measurement, i.e., age (± 5 years), sex, duration of dialysis (± 3 months), diabetes mellitus, cardiovascular disease, stroke, chronic liver disease and propensity score for PD (± 0.05), of which the independent variables included hypertension, hyperlipidemia, and rheumatological disease (Fig. 1). Similarly, these characteristics were chosen because of the association with QOL in dialysis patients [47][48][49][50] . In addition, a wider range of age-matching criteria was applied in the cross-sectional samples because of the relatively homogenous distribution of utility values within each five-year strata (Supplementary Table 2).
Estimation of LE of dialysis population. After the 1:1 matching from the national dialysis cohort, a total of 8,570 incident HD and PD patients were identified. We first estimated the survival function of these subjects using the Kaplan-Meier method during the follow-up period of 14 years. The lifetime survival function was then extrapolated up to 720 months by a semiparametric method under the assumption of constant excess hazard [51][52][53] from the end of the 14 th year. The bootstrap method of utilizing 100 repeated samples was performed to obtain the standard error of means. The estimation was achieved using the iSQoL software 54 . The accuracy and application of this method have been validated in various studies in cohorts of different illnesses [55][56][57][58] , including dialysis populations 59-61 . Estimation of QALE of the matched pairs from the national dialysis cohort. We multiplied the survival probability at each time point (dialysis vintage) with the mean utility value of QOL and summed up throughout lifetime estimates to obtain the QALE (Figs 1 and 2). Kernel smoothing using a nearest neighbor approach was applied to average the nearest 10% of the utility values of the interviewed dialysis patients. The sample size should be at least over 50 for a random sample 62 . After this, the lifetime survival function was adjusted by the corresponding mean QOL function to obtain a quality-adjusted survival curve, of which the sum of the area under this curve was the QALE of dialysis patients.

Estimation of lifetime healthcare expenditures, cost per QALY and incremental costeffectiveness ratio (ICER) for HD and PD patients.
We retrieved the total monthly medical costs reimbursed by the NHI, which were available during the period of 1998-2010 and included all healthcare expenditures for the inpatient and outpatient services, medications, examinations, and intervention procedures (e.g., dialysis), for every enrolled incident dialysis subject (n = 8,570) throughout this time period. That is, all medical costs after the initiation of dialysis therapy (including those for treating comorbidities, complications, and so on) were included, with the exception of the costs of transportation and hiring special caregivers. Summarizing any cost of medical service after the initiation of dialysis therapy can help us to clarify which dialysis modality induces higher subsequent costs, including the dialysis modalities themselves and their associated complications. These monetary values were first adjusted to those for 2010 according to the Consumer Price Index, and the average monthly cost was calculated for HD or PD for each month after beginning dialysis. The extrapolated NHI expenditures after 2010 were annually discounted at 3%, as recommended by the World Health Organization's project entitled Choosing Interventions that are Cost-Effective (WHO-CHOICE) 63 . The total average monthly expenditures were then multiplied by the monthly survival probabilities for HD or PD and summed up throughout lifetime to estimate the lifetime healthcare expenditures. The equation for ICER (PD to HD) was calculated as follows: [Lifetime cost for (PD -HD)]/[QALE of (PD -HD)].

Statistical analysis.
Continuous variables were summarized as means ± standard deviations. Comparisons of continuous variables were assessed by using the Student's t test and/or the Kruskal-Wallis test. The categorical variables were presented as the number of cases (and percentages), and the comparison of differences or trends between various groups was carried out using the Chi-square test, Fisher's exact test or Mental-Haenszel Chi-square test. All statistical analyses were performed with SAS version 9.4 (SAS Institute, Cary, NC, USA). A two-sided p value < 0.05 was considered statistically significant.

Results
Baseline characteristics of the patients from the national dialysis cohort and the cross-sectional samples. We identified 71,796 incident cases of ESRD during 1998-2010 after exclusion. After matching with individual characteristics and propensity scores, we recruited 4,285 matched-pairs of incident HD and PD patients from the national dialysis cohort (Fig. 1). Among 2,184 prevalent patients non-selectively invited from 12 dialysis centers, 295 patients refused the EQ-5D survey. After excluding patients aged under 18 (n = 28), having dialysis mode switching or renal transplantation (n = 28) and having malignancy (n = 146), 1,403 HD and 284 PD patients were entered into the final analysis. Table 1 summarizes the differences in baseline characteristics among the national dialysis cohort, responders and non-responders of the cross-sectional samples. The responders from the cross-sectional samples were more likely to be younger and with fewer comorbidities than those of the national dialysis cohort in either HD or PD groups. Compared with the non-responders, the responders of HD patients were more likely to have lower proportions of education levels at "did not finish school" or "get married" (66.6% versus 77.1%, P < 0.001), having a lower "duration on dialysis" (74.3 versus 89.5 months, P < 0.001), and to show higher prevalences of hypertension (61.5% versus 52.1%, P = 0.003) and hyperlipidemia (6.2% versus 2.1%, P = 0.005). In the PD patient group, the distributions of baseline characteristics between responders and non-responders were similar.
Baseline characteristics of the patients from the national dialysis cohort and the cross-sectional samples after matching for clinical characteristics and propensity scores. Table 2 showed there were no statistically significant differences between the matched pairs of incident HD and PD patients selected from the national dialysis cohort in age, sex, dialysis duration, urbanization and most of the comorbidities. However, PD patients had a lower prevalence of peripheral vascular diseases (9.4% versus 12.7%, standardized differences [di]:10.54), but showed slightly higher prevalences of hypertension (89.5% versus 85.8%, di:11.26) and hyperlipidemia (53.7% versus 46.7%, di:14.03). Among the matched pairs of prevalent dialysis individuals from the cross-sectional samples, there were no statistically significant differences in most of the clinical characteristics, except in dialysis duration, age, education level and employment status.

Distribution of EQ-5D scores for each domain, utility values and EQ-5D visual analogue scales in prevalent ESRD patients on HD and PD therapy.
In the EQ-5D survey among the total cross-sectional samples of prevalent dialysis patients, most PD patients reported no problems in mobility (86.6%), self-care (90.9%), usual activities (83.5%), pain/discomfort (68.0%) or anxiety/depression (75.0%) ( in mobility, usual activity and pain/discomfort. The mean EQ-5D values calculated by an equation established by Taiwan (TW) 64 were also higher in patients with PD than with HD (0.90 versus 0.81, respectively). To facilitate an international comparison, a sensitivity analysis by replacing the utility values estimated by the United Kingdom (UK)-based equation 43 and the scores of the visual analogue scale (0-100) showed the same trend. After matching, HD patients had slightly higher proportions reporting no problems in all domains than PD patients, except in the domain of anxiety/depression. When converting the utility values by the TW-based equation, the mean utilities were nearly equal for both HD and PD patients (0.88 versus 0.89, respectively). Applying the utility values estimated by the UK-based equation and the scores of the visual analogue scale also showed no differences between the matched HD and PD patients. In addition, we added employment status as one more criterion for matching and re-ran the analysis, which resulted in 144 pairs for HD and PD with utility values based on the TW function of 0.88 and 0.88, respectively (Supplementary Tables 3 and 4 for about 90% of the total cost, and PD was almost always less costly with regard to monthly total healthcare expenditures, as shown in Fig. 3. This indicates a higher dialysis cost and/or higher frequency of outpatient visits for HD, and the complications requiring hospitalization in PD were not more severe than those in HD throughout the observation period. After adjusting for a 3% annual discount rate for both denominator and numerator, the cost-effectiveness ratios for PD and HD therapy were 13,681 and 16,643 USD per QALY, indicating savings of nearly 3,000 USD per QALY for PD compared with HD. The cost per QALY of both HD and PD represented 0.74and 0.90-fold the value of gross domestic product per capita of Taiwan in 2010, which indicates dialysis therapy is cost-effective, in line with the suggestion of the WHO-CHOICE project 63 . The estimated ICER of PD to HD was − 50,858, which further supports the dominance of PD in the cost-effectiveness analysis.
Sensitivity Analysis. The sensitivity analysis, performed using 2:1 matching, showed largely the same results as those observed in the 1:1 matching. The distributions of baseline characteristics between HD and PD after matching were almost the same (Supplementary Table 5). Although there was no significant difference for survival between HD and PD during the follow-up period (Supplementary Figure 1), the estimated LE was longer for HD than for PD (

Discussion
This study adopted a matching design for a nationwide incident dialysis population followed for 14 years and comprehensively estimated lifetime survival function, costs, and dynamic changes of QOL, which might solve most of the potential limitations of previous studies. For patients without switching modes, our results show that PD is more cost-effective than HD if the dialysis cost of HD is higher than that of PD. This is strongly corroborated by the following arguments: First, since we closely matched PD and HD patients on sex, age, major    comorbidities and propensity scores, this study controlled as many potential confounders as possible that might account for the differences of survival and QOL functions and preference for dialysis modalities, which may have improved the internal validity. Second, the application of utility values throughout duration-to-date with a kernel smoother could account for the dynamic change in QOL over time. Third, since all dialysis patients registered in the catastrophic illnesses in Taiwan can be waived from co-payment, and the total follow-up period is more than a decade, the total healthcare expenditures estimated in this study would be very comprehensive and information bias would be minimized for both PD and HD. Since the resulting survival and QOL functions for PD and HD were almost the same after matching (Table 4), the major differences in the cost-effectiveness ratios mainly came from the different costs of different modalities and the associated complications. As the dialysis cost under PD is less expensive than HD in Taiwan, PD is dominant over HD with regard to cost-effectiveness, as is also seen in Fig. 3 (and Supplementary Figure 2). Numerous studies have been conducted to compare the survival outcomes between patients receiving HD and PD [14][15][16][17][18][19][20][21][22][23][24][25][26][27] . Although the messages from these studies were mixed and sometimes confusing, the majority indicated that HD and PD conferred similar survival rates, despite some differences that existed within the distinct subgroups of the patients 65 . In our study, the comparison of the survival function of the 1:1 matched pairs from the national dialysis cohort during the 14-year follow-up period revealed no statistical difference, as evaluated by the Kaplan-Meier method and log-rank test (Supplementary Figure 3). The sensitivity analysis of applying 1:2 matching showed consistent results (Supplementary Figure 1). Therefore, our study corroborates the current consensus opinion that PD was associated with a similar risk of mortality compared to HD.
The underutilization of PD therapy has been noted for decades. One of the possible reasons for this might be related to the debates with regard to long-term survival and secular change in QOL after initiation of HD or PD. In addition, the costs for PD are generally lower than HD in most countries 2 . However, the possible selection bias of HD and PD patients may account for this difference. In contrast to previous studies, we have taken advantage of the inter-linkages among different national databases and applied close matching as a major strategy to control the above biases. Given that we have shown the equivalence of survival and QOL after the close matching of HD with PD and the less costly nature of PD (Table 4, Fig. 3 and Supplementary Figure 2), we conclude that PD is dominant in the cost-effectiveness analysis. As 71-76% of dialysis patients are suitable for PD therapy [66][67][68] , our evidence can be used for policy makers to reconstruct the imbursement policy for the promotion of PD in clinical practice for those patients without contraindication.
The heavy financial burden of dialysis therapy on health care systems has led to numerous economic evaluations of different dialysis modalities 7,13,69,70 , and a Markov model is commonly used in evaluating the cost-effectiveness of these approaches. One of the basic assumptions of Markov model is the use of fixed sets of transition probabilities, health care costs, and QOL, which are usually derived from the data of a short follow-up period. In reality, these parameters are time-dependent [71][72][73][74] , which would increase the inaccuracy of the resulting estimations, and a sensitivity analysis is always needed to cover the uncertainty. As there have been no randomized control trials with a large sample for HD and PD 29 , the parameters applied in the Markov model are generally subject to the selection bias of HD and PD in the general dialysis population. In this study, we directly measured the monthly survival rates and average costs of medical services after initiation of dialysis for more than a decade ( Fig. 3 and Supplementary Figure 3). In addition, the kernel-smoothing average estimation method takes the dynamic change in QOL over time into account 62 (Supplementary Figure 4). As such, our estimates might more accurately catch the real patient-centered conditions and be closer to the actual situation.
This study showed a generally higher mean EQ-5D utility value of both prevalent HD and PD (0.83 and 0.90, respectively) than those of the previously published reports (HD group: 0.44− 0.71; PD group: 0.53-0.72) 35,75 . Such a difference might be related to the younger age and fewer comorbidities in our cross-sectional sample (Table 1), and the fact that there is no copayment for dialysis under the NHI system of Taiwan, which may reduce the financial pressure and psychological stress of dialysis patients. In fact, one previous study from Taiwan revealed that the scores for the item "health and social care availability" and the environment domain of WHOQOL-BREF were even better in HD patients than in those of healthy referents 50 , which indicates the satisfaction of dialysis patients with the NHI system. Nonetheless, a Taiwanese nationwide EQ-5D survey was performed in 2009, which showed a significantly higher average utility value of the general population of 0.931 after matching sex and age with our dialysis sample, or, a mean utility of 0.824 (Supplementary Table 7). Furthermore, we conducted this survey during 2011-2, and so the advances in dialysis technology and medical care over time and different cultural systems 76 may also be responsible for the differences. Since both national PD and HD cohorts were generally older and showed higher proportions of comorbidities than those of the cross-sectional samples (Tables 1 and 2), the utility values in this study may be overestimated, and one must be cautious in making generalizations from our study results to the whole dialysis population.
Because of low copayment, convenience and nearly full coverage of all medical services, most citizens in Taiwan are strongly adherent to the NHI program, and it provides the government and our team an opportunity to comprehensively estimate the total healthcare expenditures based on the reimbursement database. The 2012 total healthcare expenditure per capita in Taiwan was 1,350 USD 77 , which is only 15.3% of the value in the United States, 28.6% of Germany, 28.2% of Japan and 37.6% of the United Kingdom 78 . As we have also taken the 3% discount rate into consideration, the final estimated lifetime healthcare expenditure and cost-per-QALY in our study appears much lower than those seen in these high income countries. Since the average cost for each PD service is usually lower than that of HD in many high income countries, we would expect that the resulting difference in lifetime healthcare expenditures and costs-per-QALY between PD and HD would also be higher in these countries.
There are several limitations to our analyses, as follows. First, we excluded patients who received transplantation or shifting dialysis and performed individual matching to control differences in risk-factor distributions of survival function and QOL across different dialysis modalities. Though this process may improve the internal Scientific RepoRts | 6:30266 | DOI: 10.1038/srep30266 validity of the comparison between PD and HD, it might make it less applicable to generalize to the whole dialysis population; however, the detailed data in Tables 1 and 2 provide a basis for valid comparison and possible adjustment. Since the above procedures only excluded 6.6% and 1.3% of the national cohort with ESRD and patients eligible for assessing QOL, respectively, they would at most lead to a minimal bias and would not change the conclusion of PD dominance. Second, compared with the 1:1 matched national cohorts (n = 4285 in Table 2), the matched cross sectional sample (n = 179) were slightly younger in age and had lower prevalence rates of comorbidities (including diabetes, cardiovascular diseases, and liver cirrhosis), which could possibly lead to the overestimation of the QOL values for both PD and HD. Since these values were comparable for both dialyses, they would not likely change the conclusion of PD dominance. In contrast, our sensitivity analysis of adding employment status as one more criterion for matching resulted in a small sample size (n = 144), which also showed similar results (or, 0.88 vs. 0.88 for HD vs. PD, as in Supplementary Tables 3 and 4). Therefore, we expect that any change in utility values resulting from the selection criteria that improved the validity of the comparison would be unlikely to affect the conclusion of PD dominance. Further, we assumed that the values of QOL were constant and the same as those at the end of the follow-up period in order to estimate the value of QOL measured beyond the longest follow-up period in our cross-sectional samples. This could lead to the overestimation of QOL within this period because of the gradual decline of QOL during the aging process. However, the estimation bias would be minimized with the concomitant decrease in the survival rate after 14 years of dialysis. Third, the estimated expenditures did not include costs of family care and opportunity costs of the patients. In light of the flexible dialysate exchange program, patients on PD are more likely to keep their employment status and spend less time and money on transportation to dialysis centers. Therefore, the difference in costs between HD and PD would be larger if we took the opportunity costs into account. Fourth, although we have tried our best to match major comorbidities and propensity scores in this study, we still cannot completely rule out residual confounding resulted from choices and decisions involved in the selection process of dialysis modalities. For example, we do not have data on body mass index or hemoglobin level available in the two databases for our study, although both factors are related to survival function, healthcare expenditure and QOL. Besides, the lack of information with regard to decision making for dialysis modalities, such as the incentives offered by the reimbursement policy or personal lifestyle factors, might also leave some residual confounding and/or selection bias.
Despite its limitations, this study shows that the survival and QOL functions were almost the same for HD and PD after matching. The superiority of PD in this study resulted from its generally lower overall cost compared with HD. Since the total healthcare expenditures, including costs of inpatient and outpatient services, for dialysis modalities are mainly driven by the costs of the dialysis procedure itself 2 , our conclusions can probably be generalized to countries with higher expenditures for HD than PD per dialysis session. The potential saving of opportunity costs for PD further increases its superiority over HD, especially among those patients with regular employment. However, future studies are warranted to explore the cost-effectiveness of shifting dialysis modalities.