We conducted a cost-utility analysis using individual patient data to compare the standard alternative (observed) with the alternative of a new treatment (simulated) by estimating quality-adjusted life years (QALY) and costs.9 This fictitious economic evaluation compared the standard- of-care (Anatomical Therapeutic Chemical [ATC] ATC code XXX1), with a new treatment (ATC code XXX2) which has not yet been used in clinical practice for individuals with conditions that increase cardiovascular risk. Our health technology assessment scenario applies once the new medicine has been authorised, but its implementation in clinical practice is still pending. We derived treatment-related risk reductions for stroke, acute myocardial infarction (AMI), and death from clinical trial in terms of hazard ratios (HR) for each event. Since our follow-up period ends on 31 March 2023, any events occurring after this date will be considered censored on the right, as the observed time is less than the event time. We estimated extrapolated QALY and costs by assuming that the subsequent trajectory will be the same as that of the equivalent population, adjusted for excess risks after stroke and AMI and sequelae due to stroke. Figure 1 show the computational workflow for estimating individual costs and QALY according to the timing of events (stroke, AMI and death).

Figure 1.

Flow diagram of observed population, showing events and calculations of costs and quality-adjusted life years. AMI: acute myocardial infarction; HR: hazard ratio; SES: socioeconomic status; QALY: quality-adjusted life years.

We adopted a lifetime horizon. A 3% discount rate was applied to calculate the incremental cost and effectiveness of the treatment. The database contains information on whether each individual has experienced an event, as well as the date on which it occurred. As the results of the existing treatment are already available, we can calculate its QALY and costs directly. In contrast, we need to create a simulated alternative with results corresponding to the new treatment. The effect of the new treatment has only been taken into account during the follow-up period. Consequently, the cost of treatment is only allocated during this period. The extrapolated life expectancy is adjusted according to events that lead to death (survival 0) and those involving the costs of dependency and disutilities.

The results section provides an overview of the files in Excel and R, which are included in the Supplementary material, to help in understanding the exercise.

The Supplementary material contains additional methodological elements, such as a description of advanced methods not covered in this tutorial (SM1), probabilistic sensitivity analysis (SM2), rules for managing competitive risks in the allocation of events and survival times (SM3), elements of internal and external model validation (SM4), and the process for selecting the best model in a parametric survival analysis (SM5).

The dataset included individuals who were diagnosed with conditions that increase cardiovascular risk and who were treated with an ATC code XXX1 drug between 1 January 2020 and 30 June 2021. These individuals were identified by their first prescription with this code. The dataset was obtained from a real-data study, but was then modified and anonymised to enable us to maintain the correlation between variables for each individual. The follow-up period for collecting data on cardiovascular events and deaths ended on 31 March 2023. These conditions included high blood pressure, high cholesterol and atrial fibrillation. Their characteristics according to sex, socioeconomic status (SES) and the Charlson comorbidity index (CCI)10 are shown in Table 1. SES was also included to analyse potential inequities. Life expectancy data are provided to estimate the long-term impact of lives saved on effectiveness. The life expectancy of survivors is assumed to be the same as that of people of the same age and sex with the excess risk of death associated to stroke and AMI.

The societal perspective was employed by including both healthcare and social costs related to dependency. The standard treatment (ATC XXX1) costs €1,200 per year. New drug treatment (ATC XXX2) costs €2,500 per year. Both treatments were administered until the end of the follow-up period. Hospital costs for stroke vary according to length of stay, averaging €5,017, with a standard deviation of €2,050. For AMI, the average is €6,183, with a standard deviation of €2,364. Fifty percent of stroke patients have disabling conditions requiring treatment, resulting in formal annual social costs of €12,000 until death. Given that the level of autonomy of stroke patients is not recorded, determining which patients will experience sequelae is not possible. Therefore, we allocate €6,000 per year to all stroke cases.

As the treatment in question is not yet in clinical use, we drew information regarding its effectiveness from published clinical trials. For the purposes of our analysis, we will assume that these trials indicate an HR of 0.8 for all three events. In terms of survival, stroke and AMI events carry a competing risk of death. The key intermediate measure required for estimating QALY is the time taken to reach the three events (stroke, AMI and death) with the new and standard treatment alternatives. The times for the standard alternative are already known from the observed dataset. The challenge lies in estimating the times in the new alternative via a simulated dataset. To achieve this goal, it is necessary to estimate the times to three events —stroke, AMI and death— in the alternative scenario in which the new treatment has been implemented.

Given that the risk of complications and death depends on age, the impact of the two treatment alternatives must be assessed for each patient individually and extrapolated over their life expectancy. Quality of life in terms of utility values is disaggregated by decade, SES and sex and taken from the EuroQol-5D-3 L included in the 2012 National Health Survey (Table in Excel file).11–13 The patients’ disutility was 0.1786 after a stroke and 0.0365 after AMI.11

The remaining life expectancy was calculated for each individual with one parametric survival function for men and another for women based on year-by-year mortality rates obtained, from the National Statistics Institute (INE) tables for 2019 (Table SM1 in Supplementary material) following the procedure described by Arrospide et al.14. The life expectancies were validated with the expected survival values estimated by the INE in 2019 (see Table SM5 in Supplementary material). Life expectancy was adjusted for each individual by including a HR if they had suffered a stroke (men=1.6; women=1.2)15 or an AMI (men=2.0; women=2.0).16 Individuals who died had zero life expectancy. To avoid stochastic uncertainty, the same random parameter included in the function was assigned to each individual in both alternatives. This approach explicitly adjusts for comorbidities by integrating excess mortality risk into the survival function and was previously used in a cost-effectiveness analysis of a COVID-19 vaccination programme.17

This study was performed in line with the principles of the Declaration of Helsinki and received approval on May 12, 2021 from the Ethics and Clinical Research Committee of Euskadi (study code PI2021085) which waived the informed consent as all data were anonymized.

The exercise files in Excel and R (see Supplementary Material) are available in the OSF repository (https://osf.io/gyuq9/overview?view_only=748bc600f59646f1b144e21458c50a42). The description of the files can be found in Table 2.

The electronic medical records for the 1710 participants provided the following variables: treatment group: standard treatment; age: age; sex; SES: low/high as defined by the copayment category); work activity: pensioner/worker; Charlson_factor: CCI; stroke event: yes/no; AMI event: yes/no; death event: yes/no; cost of hospitalization from analytical accounting (variable Coste_proc); time to stroke event (Follow-up_stroke): time from start of treatment to stroke event or end of follow-up if no event; time to AMI event (Follow-up_AMI): time from start of treatment to AMI event or end of follow-up if no event; time to death event (Follow-up_death): time from start of treatment to death event or end of follow-up if no event; time to total follow-up (Total_follow-up): time from the start of treatment to the end of the study follow-up. If the patient has no events, these four times are equal. The time until each event is shown in days.

We derived utilities disaggregated by sex and SES from the 2012 Spanish Health Survey.11 The table shows the utilities per decade and, for each life expectancy, calculates the average utility according to how long a person will remain in each decade.

The follow-up period ends on 31 March 2023 for patients who are still alive, or on the date of their death. The estimated times to event, as calculated by the parametric functions, were converted into final times to event in the simulated dataset according to the established rules for addressing competing risks. These rules are fully described in the Supplementary Material. The simulated time to death is the estimated time if it is shorter than the final observed time for each individual. Conversely, the simulated time is the same as the observed time. For both stroke and AMI, an event is deemed to have occurred if the estimated time is shorter than the end of the follow-up period. This function also enabled the HR of each event with the new treatment (HR=0.8) to be included. The ‘Dataset_survival_times_to_event_simulated’ file displays the survival analyses and estimated event times. Assuming the same HR for all three events is unrealistic, but this simplified exercise is intended to facilitate reproduction.

The only change made to the observed dataset to create the simulated dataset was the substitution of events and times until the event, achieved by applying the rules to address the competing nature of the values estimated by the parametric functions. Therefore, as they depend on events, the QALY and costs also change. The cost of treatment also changes to account for the new drug.

First, we analysed the time to each event for the population treated with the conventional drug using parametric survival models including as covariates age, sex and CCI. This enables us to select three age-adjusted survival functions depending on the lowest Akaike information criterion (AIC) value. Second, we apply inverse Gompertz (death) and Weibull (stroke and AMI) functions to estimate the adjusted time to event in the simulated dataset. To achieve this, we solve the inverse function of the standard probability survival function with time to event as the outcome variable. Caro et al.6 presented formulas for these ‘quantile functions’ on page 115 of their book. A random parameter is added to incorporate first-order or stochastic uncertainty.18 The second step is to calculate the times to event for the simulated alternative, incorporating the HR into the function.

The formulas are provided in Excel and R files. Caro et al.6 described the combination of inputs and equations from different sources as incorporating external data on intervention effects.

Life expectancy is assigned according to age, sex, and CCI and calculated in the file “Validacion supervivencia Gompertz” (using the Life_expectancy variable).

Life expectancy is adjusted by discounting via the continuous formula (Disc_Life_Expectancy variable). Cases of death have a value of 0 in this variable and are assigned by filtering.

To assign utilities to individuals according to their age, we must consider that the utility differs for each decade. According to age, life expectancy and the utilities for each decade according to sex, the average utility by age and sex is calculated in the life expectancy sheet (Utility_Table).

To make the calculations easier to understand, the following distinction is made:

QALYs_disc_Extrap: post-follow-up or extrapolated QALYs. Discount-adjusted QALY (the ‘QALY_disc’ variable) are calculated by multiplying the ‘Disc_Life_Spe_disc’ variable by the mean utility for each age, taken from the ‘Utility_Table’. For the deceased, a value of 0 QALYs is assigned.

QALYs_disc_follow-up is the number of QALY to be calculated during the follow-up of the study. This accounts for the period before and after the stroke, as disutility needs to be applied after the stroke. These QALYs are not discounted, as the periods have already occurred19.

Finally, both QALY are summed to yield the total QALY (QALY_total).

We distinguish between three types of cost:

• •

  Treatment costs (Cost_tratam): these costs are calculated by multiplying the adjusted life expectancy+the duration of follow-up during the study by the cost of treatment according to group.

• •

  Costs of stroke sequelae (Cost_dependency): the poststroke follow-up duration and the discounted life expectancy duration are added according to the continuous formula (variable Disc_Espe_life) and multiplied by the annual sequelae cost (C_stroke_sequelae). A cost of 0 is assigned to deceased patients.

• •

  Hospital cost of events (Coste_hospital_eventos): add the costs of stroke and AMI with no discount.

Finally, the three costs are added together (variable Total_cost_disc).

Once the QALY and costs have been calculated in the observed dataset, the dataset is modified to estimate the QALYs and costs in the simulated dataset. To achieve this, the times to event are simply changed according to the previously described procedure. Once the times have changed, the event assignments are recalculated. Each case is analysed separately according to the rules described to deal the competing risks.

The dataset is analysed by calculating the incremental cost-effectiveness ratio (ICER) according to group, sex, and SES (Table 3). To achieve this, the two datasets (observed and simulated) are combined in a single file by adding the treatment variable. The results are shown in the Statistical Analysis Solved Exercise 3 file, which contains the ICER calculation disaggregated by sex and SES, as well as the cost-effectiveness plane. These quantitative results are purely fictitious and are provided to demonstrate how to complete the exercise using both Excel and R.

Further statistical analysis can be performed by adjusting for age, sex, SES and Charlson variables via SUR20 or separate GLM21 regression for the numerator and denominator. Glick's book21 serves as a reference, as his approach to economic evaluation in clinical trials can be applied to the resulting database. We applied bootstrapping to assess the uncertainty related to the sample.22 The cost-effectiveness plane generated through bootstrapping in the file Ejercicio_3_solucion_25_Excel_nuevo is shown in Figure 2.

Figure 2.

Cost-effectiveness plane obtained through bootstrapping.