### Intervention

The RHP intervention is based on case-management of CKD patients through the frequent and regular observance of the eGFR and microalbuminuria levels to prevent the progression of the renal disease, control of comorbidities status (diabetes and hypertension), and promotion of healthy lifestyles (better nutrition and exercise habits) [7]. The RHP defined the frequency of office appointments, testing, and rules for referral to specialized care. Conversely, patients in the standard of care did not have frequent follow up. They obtained outpatient visits on-demand and subject to availability. Additional details about the characteristics of care in the intervention have been provided in previous published studies [7], as well as the main differences between the intervention and control groups [8].

### Setting

The RHP was implemented in the Hospital E. Rebagliati Network. This includes the national Hospital E. Rebagliati, that receives patients from all over the country once their case is complex enough and warrants a referral. In addition, it includes primary care facilities with local reach [11].

### Target population

Target patients were adults 18 years old or above, with a diagnosis of CKD. Patients were classified in 5 states of disease progression according to the KDIGO guidelines [12] which determines specific aspects of the intervention delivery. In each visit, they received attention from the physician, nurse, and nutritionist, and performed an eGFR test for microalbuminuria detection. Additionally, some patients received social assistant and psychology care. Patients were required to visit a primary care facility in a variant frequency depending on their CKD stage: once a year for stages 1 and 2, twice for stages 3a and b, and three times for patients in stage 4. Patients in the early stages of the CKD (1–3a) must receive attention in primary care, while late-stage patients receive attention in specialized facilities. The analysis was performed for the entire cohort, no sub-groups were defined.

### Modelling approach

We performed a deterministic CEA comparing the cost and health consequences of the RHP to the standard of care. We used a Markov model defined by three health states: CKD (all stages), dialysis and death (Fig. 1). This is a compartmental model where the cohort transits across the health states at rates defined by transition probabilities; the probability of moving across states is given only by the current state, and hence no historic information is used [13].

For consistency with the time-to-event outcomes, the transition probabilities in the Markov model were estimated from fitted survival curves and informed by the previous epidemiological study that assessed the effectiveness of the RHP. We found two parametric survival curves to define the transitions from the “CKD” compartment to “Dialysis” and “Death” in the control group. For each survival curve, we identified a set of potential distributions including Weibull, exponential, logistic, log-normal, and log-logistic. The best fit to the observed data was the distribution with the lowest Akaike and Bayesian Information Criteria (AIC and BIC, respectively). In addition, we used a visual representation of the data to assess how well the curves fi the observed data.

Using the results of the selected distribution we estimate the lambda and gamma parameters as: \(1/\left[ {\exp \left( {intercept} \right) \wedge \left( {1/scale} \right)} \right]\) and \(1/scale\), respectively. These two parameters allowed us to estimate the survival probabilities from CKD to dialysis and death for each cycle in the usual care following this formula: \(1 - \exp \left[ {lambda \times \left( {cycle - 1} \right)^{gamma} - cycle^{gamma} } \right].\) Thus, this is the probability of transition from the first health state to the other ones, expressed as a function of the cycle that accounts for the increasing risk to event over time. The time scale in the previous study was days, in this study we used years for convenience given the time-horizon of the analysis. We fitted survival curves under both scales obtaining very similar results. We found that the Weibull distribution had the best fit among the evaluated parametric curves for both progressions to dialysis and mortality. For event-free survival the estimated lambda and gamma were 0.008 and 0.938, respectively; and 0.043 and 1.143, for the mortality curve. The numerical and graphical results of the calibration process for both curves are display in the Additional file 1.

For the intervention, we adjusted the survival curves by the estimated treatment effects; the hazard ratios found in the epidemiological study: 58% (HR = 0.42, 95% CI 0.21, 0.71) reduction in risk to progress to dialysis and no change in mortality risk (HR = 1, 95% CI 0.88, 1.13) [8]. Given that the target population of the epidemiological study and the CEA are the same, these treatment effect estimates have enough internal validity to create reliable results. However, we given the inherit uncertainty of the estimates, we will include both treatment effects in the sensitivity analysis.

The probability of transition from dialysis to death was estimated through literature review. We found two studies that estimated mortality rates after RTT initiation. One reported a survival probability of 95% for the first year of dialysis, 91% for the second, and 88% for the third [14]. Since our model only considers people that just started in dialysis, and the model has a memoryless property, we use a fixed rate of 5% to estimate this transition.

Considering that the mean age of the population of interest is around 60 years old [8], our model ran for 30 cycles, where each cycle represents one year, to accrue the lifetime health and economic consequences of the intervention. According to country-specific life tables, after 90 years of age, less than 10% of the Peruvian population would still be alive [15]. Cost and outcomes were observed at the end of each cycle. We used a simulated cohort of 1000 people for each strategy, replicating the same CKD stage distribution and diabetes prevalence as in the original dataset for representative purposes and internal validity.

### Cost assessment

Costs were estimated from the payer perspective, EsSalud, considering direct medical costs for both alternatives. We considered all costs faced by the payer to provide treatment in one year. For the intervention, we considered the costs of nephroprotection treatment, outpatient visits, and laboratory tests. These costs are not homogeneous across stages of CKD but vary due to frequency of provision (hospital visits and lab tests), and type of facility (primary care facilities have cheaper provision costs than more specialized ones), and patients with diabetes receive additionally glycosylated hemoglobin tests. We also included the implementation cost of the intervention, including a first investment that includes the time utilized by the Renal Health Unit to develop the intervention, protocol, guidelines, personnel training, and a yearly operational cost that includes the time spent by the RHP team identifying, testing, referring, and keeping accurate records of the patients.

While the patients undergoing intervention followed a treatment protocol and were closely followed-up, patients on the standard of care receive on-demand outpatient visits. Testing was subject to medical indication. There was not a fixed frequency for neither of those services, contrary to the intervention. Thus, healthcare utilization would depend on the patients’ behavior and availability. Given the inherent randomness of this healthcare utilization variables, we decided to base our estimates of frequency of care reception on expert consultation. The nephrologists from the Renal Health Unit provided us with their best-educated guess of the number of outpatient visits that a regular CKD patient receives in one year. We also added the cost of one laboratory test per year without differentiation for diabetes condition or facility in which it would take place. This allowed us to keep the usual care costs low, to obtain conservative estimations for the incremental cost-effectiveness ratio (ICER).

The treatment costs varied across CKD stages in both alternatives, and across diabetes status in the RHP. Given that the Markov model included one compartment for all CKD stages, we used a weighted average to estimate the annual treatment cost per patient in each alternative. The weight was defined by the proportion of stages and diabetic patients in the observed data. Finally, we included the annual cost of hemodialysis as the product of the cost per session plus the drugs prescribed in each one and the number of sessions in a year.

We used two sources of costing data: the EsSalud General Management Office cost report (2018), and the report of resources use specifically for the intervention from the Renal Health Unit (2014). The first one provided the unit cost per activity, while the latter gives us the number of units per activity consumed to follow and treat a regular patient in each CKD stage. We used institutional costs to reduce uncertainty around the final estimations. Data collection was conducted in local currency, Peruvian Soles (PEN), while results are presented in United States Dollars (USD, $).

### Health consequences

We sought to compare the differences in health outcomes between the alternatives using years lived free of dialysis (YL) and Quality Adjusted Life Years (QALY), to obtain a measure of the number of person-years avoided in dialysis and the number of person-years of perfect health gained, associated with the adherence to the intervention. Form literature review we defined the utility score for CKD patients in 0.84, and for patients starting dialysis in 0.65 [16]. In the study, these scores correspond to the stage 3, the most prevalent stage in our cohort, and to the stage 5, a close approximation to patients just starting RRT.

### Analysis

We projected the costs and health outcomes of each alternative separately during 30 cycles, to capture life-time consequences. Cost and health outcomes would be discounted by an annual rate of 3% to reflect the time-preferences of the economic agents. After these calculations, we aggregated the total costs, YL and QALY from each alternative and express them in per-person units.

To determine which alternative poses the highest economic value we used the ICER calculated as the difference in cost between RHP and usual care over the differences in health outcomes. Then, we can interpret the ICER as the additional cost for the payer to avoid one person-year in dialysis and to gain one QALY. Cost estimations were made in local currency but converted to USD using a fixed exchange rate of 3.3 PEN per each USD, corresponding to the annual average for 2018 [17]. We used the cost-effectiveness threshold of 1–3 times the gross domestic product (GDP) per capita, estimated in $6571 for Peru, according to the World Bank [18].

To assess the robustness of the ICER against parameter uncertainty we performed a Probabilistic Sensitivity Analysis (PSA) based on a Monte Carlo simulation of 1000 repetitions. In each repetition, the model randomly picked a value for each varying parameter considering its distribution and range of values. Each repletion represents a unique scenario for the comparison of RHP and the standard of care. Table 1 shows a summary of the parameters used in the study and the values they took in the sensitivity analysis. The range of values for each parameter was as follows: the treatment effects would take the lower and upper values of the estimated confidence interval, the costs would vary 15%, the utilities would change 10%, and the discount rate would take 0% to reflect no discounting and 5% to reflect a scenario with higher opportunity cost. We summarized the results by descriptive statistics of the incremental cost per QALY and a figure showing the incremental costs and QALYs for each simulation. We report central tendency values for the distribution of both the incremental cost and effectiveness, including mean, standard deviation (SD), relative standard error \(\left( {RSE = SD/mean \times 100} \right)\) expressed in percentages, and range (min. and max.) of values.

The analysis and assessment of parametric distributions to fit the data was performed in R studio, the CEA and PSA were developed in TreeAge^{®} [19].