MSc - FM4 - Health economics

Applied health economic evaluations using R

Authors

Affiliation

Nathanael Lutz

Bern University of Applied Sciences

Jan Taeymans

Bern University of Applied Sciences

Published

October 8, 2026

1 Introduction

This hands-on exercise provides an opportunity to apply the theoretical concepts introduced in the health economics lectures by Prof. Dr. Taeymans and reinforce key ideas through practical application. It covers the central steps of a health economic evaluation, while making some deliberate simplifications to remain manageable within the available teaching time. Nevertheless, if you work in this field in the future, the exercise can serve as a useful starting point and a resource to which you can return.

The exercise consists of two parts:

a trial-based health economic evaluation
a Markov model

I tried my best to include in each part a short recap of the relevant theory, a physiotherapy-related example, specific tasks and solutions with explanations. Model-based health economic evaluations can become technically complex rather quickly. The Markov exercises therefore focus on the fundamental structure, calculations and interpretation of these models rather than the level of detail required for a publication-quality evaluation.

We will use R because it is free, open source, and supported by an extensive user community. As you noticed, R can be used for almost every type of quantitative data management and analysis. It also facilitates transparent and reproducible research, making it a valuable and widely accessible tool.

Some basic familiarity with R will be helpful, but you are not expected to know every command by heart. We will provide the essential code and guide you through it.

Before copying, pasting, and running a piece of code, take a moment to understand what it does. Otherwise, you may obtain the correct result without understanding how you arrived at it.

You can complete this exercise in a new R script or, if you are already familiar with them, in a Markdown or Quarto document. A Markdown or Quarto document may help you keep your code, results and interpretations together and provide a clearer overview of the complete exercise.

We expect the complete set of exercises to require at least eight hours of work. Please use the scheduled on-site sessions to ask questions, discuss difficulties and obtain support.

Content

Trial based HEE
- Intervention costs
- Direct and indirect costs
- Discounting
- QALYs
- ICER, ICUR
- Cost-effectiveness plane, cost-effectiveness acceptability curve
- Uncertainty and sensitivity analysis
Markov model
- define health states
- specify transition probabilities
- construct a transition matrix
- follow a cohort over several cycles
- assign costs and utilities
- calculate total discounted costs and QALYs
- compare two strategies
References and further reading

2 Trial-based health economic evaluation

In this exercise, we evaluate a structured physiotherapy programme for people with knee osteoarthritis over 24 months. The programme combines patient education, supervised exercises and an individual home-exercise programme. The treating physiotherapists determine the number of supervised sessions according to each participant’s clinical needs and progress.

Participants in the control group receive usual care. They do not participate in the structured programme but may use healthcare services available in routine practice, including additional physiotherapy, physician consultations, medication and hospital care.

The effectiveness of the programme is measured using the KOOS Activities of Daily Living subscale (KOOS-ADL). Scores range from 0 to 100, with higher scores indicating better function. EQ-5D-5L responses are used to estimate health-state utilities and calculate quality-adjusted life years (QALYs).

The analysis is conducted from both a healthcare and a societal perspective. The societal perspective additionally includes productivity losses from paid and unpaid work.

Participants entered the study in different calendar years. Because the dataset records resource quantities rather than historical payments, we value all resources using the same 2026 unit costs. The resulting costs are therefore already expressed in a common price year. Costs and outcomes occurring after the first year are discounted.

For simplicity, direct medical resource use is limited to physiotherapy, physician visits, medication and hospitalizations. Real-world evaluations usually include more cost categories.

The dataset is artificial and was created for teaching purposes. Real trial data are often incomplete, inconsistent and affected by missing observations. Addressing these issues would add considerable complexity and distract from the main concepts covered in this exercise.

There are two data sets:

df-hee-knee-oa.rda contains the individual data of each participant. See Table 1 for details of the variables.
df-unit-costs.rda contains the unit costs for costing calculations.

You can find the data on Moodle.

Table 1: Codebook of trial based data for the health economic evaluation

Codebook for the trial dataset
Variable	Description	Type	Unit or range	Coding or definition
id	Unique participant identifier.	Character	Not applicable	One unique code per participant.
group	Randomised treatment group.	Character	Not applicable	'Control' = usual care; 'Intervention' = structured physiotherapy programme.
date_0	Date of the baseline assessment.	Date	Calendar date	Stored as an R Date
date_6	Date of the 6-month follow-up assessment.	Date	Calendar date	Stored as an R Date
date_12	Date of the 12-month follow-up assessment.	Date	Calendar date	Stored as an R Date
date_24	Date of the 24-month follow-up assessment.	Date	Calendar date	Stored as an R Date
cost_year_y1	Calendar year assigned to resource use and productivity losses occurring during the first follow-up year.	Integer	Calendar year	Applies to the period from baseline to 12 months.
cost_year_y2	Calendar year assigned to resource use and productivity losses occurring during the second follow-up year.	Integer	Calendar year	Applies to the period from 12 to 24 months.
koos_adl_0	KOOS Activities of Daily Living score at baseline.	Numeric	0–100 points	Higher values indicate better functioning in activities of daily living.
koos_adl_6	KOOS Activities of Daily Living score at 6 months.	Numeric	0–100 points	Higher values indicate better functioning in activities of daily living.
koos_adl_12	KOOS Activities of Daily Living score at 12 months.	Numeric	0–100 points	Higher values indicate better functioning in activities of daily living.
koos_adl_24	KOOS Activities of Daily Living score at 24 months.	Numeric	0–100 points	Higher values indicate better functioning in activities of daily living.
eq5d_0	EQ-5D-5L health state at baseline.	Character	Five-digit health-state code	The five digits (from 1 = no problems to 5 = extreme problems or inability) represent mobility, self-care, usual activities, pain/discomfort and anxiety/depression.
eq5d_6	EQ-5D-5L health state at 6 months.	Character	Five-digit health-state code	The five digits (from 1 = no problems to 5 = extreme problems or inability) represent mobility, self-care, usual activities, pain/discomfort and anxiety/depression.
eq5d_12	EQ-5D-5L health state at 12 months.	Character	Five-digit health-state code	The five digits (from 1 = no problems to 5 = extreme problems or inability) represent mobility, self-care, usual activities, pain/discomfort and anxiety/depression.
eq5d_24	EQ-5D-5L health state at 24 months.	Character	Five-digit health-state code	The five digits (from 1 = no problems to 5 = extreme problems or inability) represent mobility, self-care, usual activities, pain/discomfort and anxiety/depression.
intervention_sessions	Number of structured programme physiotherapy sessions attended.	Integer	Sessions	Zero for participants in the control group.
additional_physio_y1	Number of additional physiotherapy sessions received outside the study programme during Year 1.	Integer	Sessions	Covers the period from baseline to 12 months.
additional_physio_y2	Number of additional physiotherapy sessions received outside the study programme during Year 2.	Integer	Sessions	Covers the period from 12 to 24 months.
physician_visits_y1	Number of knee-related physician consultations during Year 1.	Integer	Visits	Covers the period from baseline to 12 months.
physician_visits_y2	Number of knee-related physician consultations during Year 2.	Integer	Visits	Covers the period from 12 to 24 months.
medication_months_y1	Number of months with regular use of knee-related medication during Year 1.	Integer	Months	Integer from 0 to 12.
medication_months_y2	Number of months with regular use of knee-related medication during Year 2.	Integer	Months	Integer from 0 to 12.
hospitalizations_y1	Number of knee-related hospitalizations during Year 1.	Integer	Hospital admissions	Includes knee-related inpatient treatment during the period from baseline to 12 months.
hospitalizations_y2	Number of knee-related hospitalizations during Year 2.	Integer	Hospital admissions	Includes knee-related inpatient treatment during the period from 12 to 24 months.
productivity_hours_lost_y1	Productive time completely lost because of the knee condition during Year 1.	Numeric	Hours	Includes paid and unpaid productive activities.
productivity_hours_lost_y2	Productive time completely lost because of the knee condition during Year 2.	Numeric	Hours	Includes paid and unpaid productive activities.
productivity_hours_impaired_y1	Time spent performing productive activities while impaired by the knee condition during Year 1.	Numeric	Hours	Includes paid and unpaid productive activities.
productivity_hours_impaired_y2	Time spent performing productive activities while impaired by the knee condition during Year 2.	Numeric	Hours	Includes paid and unpaid productive activities.
productivity_impairment_score_y1	Average degree of impairment while performing productive activities during Year 1.	Numeric	0–10 points	0 = no impairment; 10 = completely unable to perform the activity.
productivity_impairment_score_y2	Average degree of impairment while performing productive activities during Year 2.	Numeric	0–10 points	0 = no impairment; 10 = completely unable to perform the activity.

2.1 Import and inspect the data

Exercise

First, import the dataset and explore its contents. Set the working directory to the folder containing the data, then use load("df-hee-knee-oa.rda") to load it into R. Inspect the data using a few simple tables and plots. The functions boxplot() and table() are useful for this purpose. For example, try boxplot(koos_adl_24 ~ group, data = df.knee) or table(df.knee$group, df.knee$additional_physio_y1).

Summarise the main results from this data inspection.

There is no single correct solution to this task. Below are some examples of how to load and inspect the data.

Code

load("../data/df-hee-knee-oa.rda")

Note: the name of the data frame is df.knee.

Code

table(df.knee$group)


     Control Intervention 
         100          100

The participants are equally distributed between the two groups.

Code

table(df.knee$cost_year_y1)


2019 2020 2021 2022 2023 
  23   42   62   46   27

Participants entered the study in different calendar years. In an analysis based on historical payments, these years would be needed to adjust costs to a common price year. In this exercise, however, all resource quantities are valued directly using 2026 unit costs.

Code

table(
  df.knee$group,
  df.knee$additional_physio_y1
)

              
                0  1  2  3  4  5  6  7  8  9 10
  Control       9  9 26 18 15  8 11  2  1  0  1
  Intervention  9 20 30 21  4  9  5  1  0  1  0

The number of additional physiotherapy sessions varies between participants and groups.

Code

table(
  df.knee$group,
  df.knee$hospitalizations_y1
)

              
                0  1
  Control      91  9
  Intervention 95  5

Most participants have no hospitalization. Only a small number generate hospital costs.

Code

boxplot(
  koos_adl_0 ~ group,
  data = df.knee,
  ylab = "KOOS-ADL",
  main = "KOOS-ADL at baseline"
)

The groups have broadly similar KOOS-ADL scores at baseline.

Code

boxplot(
  koos_adl_24 ~ group,
  data = df.knee,
  ylab = "KOOS-ADL",
  main = "KOOS-ADL at 24 months"
)

At 24 months, KOOS-ADL scores tend to be higher in the intervention group.

Code

boxplot(
  medication_months_y1 ~ group,
  data = df.knee,
  ylab = "Months",
  main = "Medication use in Year 1"
)

Participants in the intervention group tend to use medication for fewer months.

Code

boxplot(
  productivity_hours_lost_y1 ~ group,
  data = df.knee,
  ylab = "Hours",
  main = "Productivity hours lost in Year 1"
)

Productivity losses are right-skewed, with a small number of participants reporting many lost hours.

Code

plot(df.knee$koos_adl_0, df.knee$koos_adl_24)

There is also a correlation between pre and post values for the KOOS.

Only by looking at the data descriptively, we expect that the intervention is more effective. Furthermore, based on the plots and tables above, we do not expect that direct medical costs and indirect costs are substantially higher in the intervention group. So from a first impression, the intervention should be cost-effective.

2.2 Intervention, direct medical and indirect costs

Healthcare utilisation can be recorded either as resource quantities or as monetary amounts. This distinction determines how costs from different calendar years are handled.

If resource quantities are available, such as the number of physiotherapy sessions, physician consultations or hospital admissions, the same set of unit costs can be applied to all participants:

\[ C_i = \sum_{j=1}^{J} q_{ij} \times p_j, \]

where:

$q_{ij}$ is the quantity of resource $j$ used by participant $i$;
$p_j$ is the unit cost of resource $j$ in the selected price year.

For example, if we apply 2026 unit costs to all observed resource quantities, every calculated cost is automatically expressed in 2026 prices. A physiotherapy session receives the same monetary value regardless of whether it occurred in 2019 or 2024. Differences in costs between participants therefore reflect differences in resource use rather than changes in prices over time.

If the dataset instead contains actual historical payments, invoices or reimbursements, these monetary amounts must first be assigned to their original price years and then adjusted to a common price year using an appropriate price index:

\[ C_i^{2026} = C_i^{t} \times \frac{\text{Price index}_{2026}} {\text{Price index}_{t}}. \]

where:

$C_i^{2026}$ is the cost for participant $i$, expressed in 2026 prices;
$C_i^{t}$ is the historical cost for participant $i$, expressed in prices from year $t$;
$\text{Price index}_{2026}$ is the value of the relevant price index in 2026;
$\text{Price index}_{t}$ is the value of the same price index in the year to which the historical cost refers;
$t$ is the original price year of the historical cost.

Why collect resource quantities?

Recording resource quantities separately from prices makes an economic evaluation more flexible. The same resource-use data can be valued using prices from different years, countries or healthcare systems without changing the underlying trial data.

In this exercise, we therefore apply one common set of unit costs to all participants. Because all unit costs refer to 2026, the calculated costs are already expressed in 2026 CHF and do not require any further adjustment for inflation.

The unit costs required for this exercise are stored in df-unit-costs.rda.

The number of consultations must be multiplied by the cost per consultation, while medication use could be valued per dose, package or treatment day. Resources should also not be counted twice. For example, physiotherapy sessions included in the intervention costs should not also be included in direct medical costs (double counting).

We calculate the following cost categories separately:

Intervention costs: resources required to provide the structured physiotherapy programme.
Direct medical costs: other healthcare resources, such as additional physiotherapy, medical consultations and medication.
Indirect costs: productivity losses resulting from absence or reduced productivity during paid or unpaid work.

Keeping these categories separate makes it possible to analyse the results from different perspectives:

The healthcare perspective includes intervention costs and direct medical costs.
The societal perspective additionally includes indirect costs.

The perspective is important because it determines which costs are relevant. An intervention may increase costs for the healthcare system but reduce productivity losses for participants and society.

There is an ongoing debate on how to value productivity. In this exercise, productivity losses are valued using the human capital approach (which is the most frequently used approach). Hours lost because of absence are multiplied by a common hourly productivity value. Reduced productivity while working is estimated using the reported hours and the participant’s impairment score (see below).

Simplifications for this exercise

Real-world costing is often considerably more detailed. Medication costs, for example, may depend on individual doses, treatment duration, package size and unused medication. Healthcare services may also have different prices depending on the provider or setting.

To keep the exercise manageable, we group direct medical resources into a small number of categories and use predefined unit costs.

Exercise

1. Load and inspect the unit costs

First, load the file containing the unit costs and inspect its contents:

load("../data/df-unit-costs.rda")

Remember that load() restores the object under the name with which it was saved. The unit-cost dataset will therefore be available as df.unit.costs - the way I like to name my R-objects 🤓

Use df.unit.costs or View(df.unit.costs) to examine the available unit costs.

For each participant, calculate the costs associated with every resource category. In general, an individual cost is calculated by multiplying the quantity of a resource used by its corresponding unit cost:

\[ C_{ij}^{2026} = q_{ij} \times p_j^{2026}, \]

where:

$C_{ij}^{2026}$ is the cost of resource $j$ for participant $i$, expressed in 2026 CHF;
$q_{ij}$ is the quantity of resource $j$ used by participant $i$;
$p_j^{2026}$ is the 2026 unit cost of resource $j$.

Because df.unit.costs contains only one row, the unit costs can be accessed directly. For example, the unit cost of one physician consultation is:

Code

df.unit.costs$physician_consultation

Physician costs during the first follow-up year can therefore be calculated as:

Code

df.knee$physician_costs_y1 <-
  df.knee$physician_visits_y1 *
  df.unit.costs$physician_consultation

2. Calculate intervention costs

Use this approach to add the following cost variables to df.knee:

intervention_costs;
additional_physio_costs_y1;
additional_physio_costs_y2;
physician_costs_y1;
physician_costs_y2;
medication_costs_y1;
medication_costs_y2;
hospitalization_costs_y1;
hospitalization_costs_y2;
absence_costs_y1;
absence_costs_y2;
presenteeism_costs_y1;
presenteeism_costs_y2.

Assume that all intervention sessions occurred during the first follow-up year. This timing assumption will become relevant when costs are discounted later in the exercise.

For absenteeism, multiply the number of productive hours completely lost by the monetary value of one productive hour:

\[ C_{\text{absence}} = h_{\text{lost}} \times v_h, \]

where:

$h_{\text{lost}}$ is the number of productive hours completely lost;
$v_h$ is the 2026 value of one productive hour.

For presenteeism, adjust the affected hours according to the reported impairment score:

\[ C_{\text{presenteeism}} = h_{\text{impaired}} \times \frac{s_{\text{impairment}}}{10} \times v_h, \]

where:

$h_{\text{impaired}}$ is the number of productive hours affected by impairment;
$s_{\text{impairment}}$ is the impairment score from 0 to 10;
$v_h$ is the 2026 value of one productive hour.

For example, an impairment score of 4 means that 40% of the affected productive time is treated as lost.

Keep all cost components separate at this stage. This will allow us to combine them later according to the healthcare and societal perspectives.

Because the same 2026 unit costs are applied to all participants, the resulting costs are already expressed in a common price year. No further adjustment for inflation is required.

1. Load and inspect the unit costs

The unit-cost file is an .rda file. The function load() restores the object under its saved name, so no assignment is required.

Code

load("../data/df-unit-costs.rda")

df.unit.costs

  price_year currency intervention_session additional_physio_session
1       2026      CHF                   55                        55
  physician_consultation medication_month hospital_admission productive_hour
1                    165               47               6800              58

The dataset contains one common set of unit costs expressed in 2026 CHF. These unit costs are applied to all participants, irrespective of the calendar year in which the resources were used.

2. Calculate intervention costs

Intervention costs are calculated by multiplying the number of intervention sessions by the cost per session.

Code

df.knee$intervention_costs <-
  df.knee$intervention_sessions *
  df.unit.costs$intervention_session

boxplot(intervention_costs ~ group, data = df.knee)

Participants in the control group have intervention costs of zero because they did not receive the structured physiotherapy programme. Intervention costs are kept separate from the costs of additional physiotherapy.

We assume that all intervention sessions occurred during the first follow-up year. This timing assumption will become relevant when costs are discounted later in the exercise.

For each year, the number of additional physiotherapy sessions is multiplied by the common 2026 cost per session.

Code

df.knee$additional_physio_costs_y1 <-
  df.knee$additional_physio_y1 *
  df.unit.costs$additional_physio_session

df.knee$additional_physio_costs_y2 <-
  df.knee$additional_physio_y2 *
  df.unit.costs$additional_physio_session

Intervention sessions must not be included in these variables. Otherwise, the same resource would be counted twice.

Physician costs are calculated by multiplying the number of consultations by the 2026 cost per consultation.

Code

df.knee$physician_costs_y1 <-
  df.knee$physician_visits_y1 *
  df.unit.costs$physician_consultation

df.knee$physician_costs_y2 <-
  df.knee$physician_visits_y2 *
  df.unit.costs$physician_consultation

Participants with more physician consultations have higher physician costs. Because the same unit cost is used for everyone, differences between participants reflect differences in the number of consultations rather than changes in prices over time.

Medication use is measured as the number of months during which a participant regularly used knee-related medication.

Code

df.knee$medication_costs_y1 <-
  df.knee$medication_months_y1 *
  df.unit.costs$medication_month

df.knee$medication_costs_y2 <-
  df.knee$medication_months_y2 *
  df.unit.costs$medication_month

This is a simplified calculation based on an average cost per month. Real medication costing would usually consider the specific medication, dose, frequency, package size and treatment duration (Jan Taeymans calls it Sherlock-Holmes-work).

Hospitalisation costs are calculated by multiplying the number of admissions by the 2026 cost per admission.

Code

df.knee$hospitalization_costs_y1 <-
  df.knee$hospitalizations_y1 *
  df.unit.costs$hospital_admission

df.knee$hospitalization_costs_y2 <-
  df.knee$hospitalizations_y2 *
  df.unit.costs$hospital_admission

hist(df.knee$hospitalization_costs_y2)

Hospitalisations are uncommon but expensive. A small number of admissions can therefore have a substantial influence on the distribution of healthcare costs (see histogram)

Absenteeism represents productive time that was completely lost because of the knee condition. The number of lost hours is multiplied by the 2026 monetary value of one productive hour.

Code

df.knee$absence_costs_y1 <-
  df.knee$productivity_hours_lost_y1 *
  df.unit.costs$productive_hour

df.knee$absence_costs_y2 <-
  df.knee$productivity_hours_lost_y2 *
  df.unit.costs$productive_hour

hist(df.knee$absence_costs_y2)

These calculations follow the human capital approach. The hourly value represents the societal value of productive time and is not the participant’s individual wage. Again, such distributions are typically right-skewed.

Presenteeism represents productive time during which a participant continued performing an activity but was impaired. Dividing the impairment score by 10 converts it into a proportion. For example, 10 impaired hours with an impairment score of 4 are treated as $10 \times 0.4 = 4$ completely lost productive hours.

Code

df.knee$presenteeism_costs_y1 <-
  df.knee$productivity_hours_impaired_y1 *
  (df.knee$productivity_impairment_score_y1 / 10) *
  df.unit.costs$productive_hour

df.knee$presenteeism_costs_y2 <-
  df.knee$productivity_hours_impaired_y2 *
  (df.knee$productivity_impairment_score_y2 / 10) *
  df.unit.costs$productive_hour

hist(df.knee$presenteeism_costs_y2)

Finally, we inspect all newly created cost variables together.

Code

cost.variables <- c(
  "intervention_costs",
  "additional_physio_costs_y1",
  "additional_physio_costs_y2",
  "physician_costs_y1",
  "physician_costs_y2",
  "medication_costs_y1",
  "medication_costs_y2",
  "hospitalization_costs_y1",
  "hospitalization_costs_y2",
  "absence_costs_y1",
  "absence_costs_y2",
  "presenteeism_costs_y1",
  "presenteeism_costs_y2"
)

summary(df.knee[cost.variables])

 intervention_costs additional_physio_costs_y1 additional_physio_costs_y2
 Min.   :  0.0      Min.   :  0.0              Min.   :  0.0             
 1st Qu.:  0.0      1st Qu.:110.0              1st Qu.: 55.0             
 Median :110.0      Median :110.0              Median :110.0             
 Mean   :258.5      Mean   :154.8              Mean   :137.8             
 3rd Qu.:495.0      3rd Qu.:220.0              3rd Qu.:165.0             
 Max.   :660.0      Max.   :550.0              Max.   :495.0             
 physician_costs_y1 physician_costs_y2 medication_costs_y1 medication_costs_y2
 Min.   :  0.0      Min.   :  0.0      Min.   :  0.0       Min.   :  0        
 1st Qu.:165.0      1st Qu.:  0.0      1st Qu.:188.0       1st Qu.:141        
 Median :165.0      Median :165.0      Median :235.0       Median :235        
 Mean   :255.8      Mean   :220.3      Mean   :251.2       Mean   :223        
 3rd Qu.:330.0      3rd Qu.:330.0      3rd Qu.:329.0       3rd Qu.:282        
 Max.   :990.0      Max.   :825.0      Max.   :564.0       Max.   :517        
 hospitalization_costs_y1 hospitalization_costs_y2 absence_costs_y1
 Min.   :   0             Min.   :    0            Min.   :   0.0  
 1st Qu.:   0             1st Qu.:    0            1st Qu.:   0.0  
 Median :   0             Median :    0            Median :   0.0  
 Mean   : 476             Mean   :  612            Mean   : 636.4  
 3rd Qu.:   0             3rd Qu.:    0            3rd Qu.: 742.4  
 Max.   :6800             Max.   :13600            Max.   :7888.0  
 absence_costs_y2  presenteeism_costs_y1 presenteeism_costs_y2
 Min.   :    0.0   Min.   :   0.0        Min.   :    0.00     
 1st Qu.:    0.0   1st Qu.:   0.0        1st Qu.:    0.00     
 Median :    0.0   Median : 332.3        Median :   95.06     
 Mean   :  672.3   Mean   : 972.3        Mean   :  943.83     
 3rd Qu.:  915.0   3rd Qu.:1473.2        3rd Qu.: 1118.72     
 Max.   :16941.8   Max.   :6559.8        Max.   :11163.49

The cost components remain separate because they will later be combined according to the analytical perspective. Because the same 2026 unit costs were applied to all participants, all calculated costs are already expressed in 2026 CHF. No further adjustment for inflation is required. We see that healthcare costs are typically right-skewed because most participants incur relatively low or moderate costs, while a small number incur very high costs.

2.3 Total costs and discounting

The individual cost components must now be combined according to the analytical perspective. Before calculating total costs, we first create annual subtotals for direct medical costs and productivity costs. Keeping the years separate allows costs occurring during the second year to be discounted.

From the healthcare perspective, we include the resources used to provide the intervention and all additional direct medical care:

\[ C_{\text{healthcare}} = C_{\text{intervention}} + C_{\text{direct medical}}. \]

From the societal perspective, we additionally include productivity losses due to absenteeism and presenteeism:

\[ C_{\text{societal}} = C_{\text{healthcare}} + C_{\text{productivity}}. \]

The choice of perspective matters because an intervention may increase costs within the healthcare system while reducing productivity losses elsewhere in society.

Costs occurring at different times are not necessarily valued equally. At baseline, costs incurred during the second year are still future costs. Discounting converts these future costs to their value at the time of the initial treatment decision:

\[ PV(C_t) = \frac{C_t}{(1+r)^t}, \]

where:

$PV(C_t)$ is the present value of the cost;
$C_t$ is the cost occurring at time $t$;
$r$ is the annual discount rate;
$t$ is the number of years between baseline and the occurrence of the cost.

In this exercise, we use an annual discount rate of 3%. We apply the following simplified timing convention:

intervention and year-1 costs are not discounted;
year-2 costs are discounted by one year.

The discounted value of year-2 costs is therefore:

\[ PV(C_{\text{year 2}}) = \frac{C_{\text{year 2}}}{1.03}. \]

For an evaluation covering several years, the discounted total cost would be:

\[ PV(C) = \sum_{t=0}^{T} \frac{C_t}{(1+r)^t}, \]

with $T$ as the final year of the evaluation.

Discounting is not cost standardisation

All costs have already been valued using 2026 unit costs and are therefore expressed at the same price level. Discounting addresses a different question: it accounts for when the costs occur relative to the treatment decision at baseline.

Exercise

The individual cost components must now be combined into annual cost categories, discounted where appropriate and then added according to the healthcare and societal perspectives.

1. Calculate annual direct medical costs

For each participant, calculate the direct medical costs separately for year 1 and year 2. Include:

additional physiotherapy;
physician consultations;
medication;
hospitalisations.

Do not include intervention costs or productivity costs in these variables.

Add the following variables to df.knee:

direct_medical_costs_y1;
direct_medical_costs_y2.

For example, direct medical costs during year 1 can be calculated as:

Code

df.knee$direct_medical_costs_y1 <-
  df.knee$additional_physio_costs_y1 +
  df.knee$physician_costs_y1 +
  df.knee$medication_costs_y1 +
  df.knee$hospitalization_costs_y1

Use the corresponding year-2 variables to calculate direct_medical_costs_y2.

2. Calculate annual productivity costs

Productivity costs consist of costs due to absenteeism and presenteeism. Calculate them separately for year 1 and year 2.

Add the following variables to df.knee:

productivity_costs_y1;
productivity_costs_y2.

The general calculation is:

\[ C_{\text{productivity}} = C_{\text{absence}} + C_{\text{presenteeism}}. \]

For example, year-1 productivity costs can be calculated as:

Code

df.knee$productivity_costs_y1 <-
  df.knee$absence_costs_y1 +
  df.knee$presenteeism_costs_y1

Use the corresponding year-2 variables to calculate productivity_costs_y2.

3. Discount year-2 costs

Use an annual discount rate of 3%. Under our simplified timing convention:

intervention and year-1 costs are not discounted;
year-2 costs are discounted by one year.

Create discounted versions of the year-2 cost categories without overwriting the original variables:

direct_medical_costs_y2_discounted;
productivity_costs_y2_discounted.

For example:

Code

discount_rate <- 0.03

df.knee$direct_medical_costs_y2_discounted <-
  df.knee$direct_medical_costs_y2 /
  (1 + discount_rate)

Use the same approach to discount productivity_costs_y2.

4. Calculate total costs from the healthcare perspective

The healthcare perspective includes:

intervention costs;
direct medical costs during year 1;
discounted direct medical costs during year 2.

Calculate the total healthcare costs for each participant and store them as:

total_costs_healthcare.

The calculation follows:

\[ C_{\text{healthcare}} = C_{\text{intervention}} + C_{\text{direct medical, year 1}} + PV(C_{\text{direct medical, year 2}}). \]

You can start the calculation as follows:

Code

df.knee$total_costs_healthcare <-
  df.knee$intervention_costs +
  # Add the year-1 and discounted year-2 direct medical costs

Make sure that intervention costs are included only once.

5. Calculate total costs from the societal perspective

The societal perspective includes all costs from the healthcare perspective and additionally includes:

productivity costs during year 1;
discounted productivity costs during year 2.

Calculate the total societal costs for each participant and store them as:

total_costs_societal.

The calculation follows:

\[ C_{\text{societal}} = C_{\text{healthcare}} + C_{\text{productivity, year 1}} + PV(C_{\text{productivity, year 2}}). \]

You can start with:

Code

df.knee$total_costs_societal <-
  df.knee$total_costs_healthcare +
  # Add the year-1 and discounted year-2 productivity costs

Using total_costs_healthcare in this calculation avoids repeating all healthcare cost components from above.

6. Inspect the total cost distributions

Create separate boxplots to compare the distributions between the intervention and control groups:

Code

boxplot(
  total_costs_healthcare ~ group,
  data = df.knee,
  ylab = "Total healthcare costs (2026 CHF)",
  xlab = "Group"
)

Adapt this code to create a corresponding boxplot for total_costs_societal.

When interpreting the plots, consider:

whether costs differ between the groups;
whether the societal perspective changes the apparent difference;
whether the distributions appear symmetric or right-skewed;
whether there are participants with particularly high costs.

There is no need to calculate incremental costs yet. At this stage, the aim is to obtain transparent participant-level total costs from both perspectives.

1. Calculate annual direct medical costs

Direct medical costs include additional physiotherapy, physician consultations, medication and hospitalisations. We calculate them separately for the two follow-up years so that year-2 costs can subsequently be discounted.

Code

df.knee$direct_medical_costs_y1 <-
  df.knee$additional_physio_costs_y1 +
  df.knee$physician_costs_y1 +
  df.knee$medication_costs_y1 +
  df.knee$hospitalization_costs_y1

df.knee$direct_medical_costs_y2 <-
  df.knee$additional_physio_costs_y2 +
  df.knee$physician_costs_y2 +
  df.knee$medication_costs_y2 +
  df.knee$hospitalization_costs_y2

The annual direct medical costs combine all healthcare services used outside the structured intervention. Hospitalisations are likely to produce some particularly high values.

2. Calculate annual productivity costs

Productivity costs combine losses due to absenteeism and presenteeism.

Code

df.knee$productivity_costs_y1 <-
  df.knee$absence_costs_y1 +
  df.knee$presenteeism_costs_y1

df.knee$productivity_costs_y2 <-
  df.knee$absence_costs_y2 +
  df.knee$presenteeism_costs_y2

These costs are relevant from the societal perspective because productive time has value regardless of whether it is lost through complete absence or reduced performance.

3. Discount year-2 costs

We use an annual discount rate of 3%. Intervention and year-1 costs remain undiscounted, while year-2 costs are divided by $1.03$.

Code

discount_rate <- 0.03

df.knee$direct_medical_costs_y2_discounted <-
  df.knee$direct_medical_costs_y2 /
  (1 + discount_rate)

df.knee$productivity_costs_y2_discounted <-
  df.knee$productivity_costs_y2 /
  (1 + discount_rate)

The discounted values are slightly lower than the original year-2 costs. With a time horizon of only two years, the difference is modest—but discounting becomes increasingly important when the time horizon grows (for example in models with a life-time horizon).

We retain the original undiscounted variables. This makes the calculation transparent and allows us to check or modify the discount rate later.

4. Calculate total costs from the healthcare perspective

The healthcare perspective includes intervention costs and direct medical costs. Productivity costs are not included.

Code

df.knee$total_costs_healthcare <-
  df.knee$intervention_costs +
  df.knee$direct_medical_costs_y1 +
  df.knee$direct_medical_costs_y2_discounted

The resulting variable contains the total healthcare cost for each participant over the complete 24-month follow-up period.

5. Calculate total costs from the societal perspective

The societal perspective includes all healthcare costs and additionally includes productivity losses.

Code

df.knee$total_costs_societal <-
  df.knee$total_costs_healthcare +
  df.knee$productivity_costs_y1 +
  df.knee$productivity_costs_y2_discounted

summary(df.knee$total_costs_societal)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  679.7  2003.1  4189.1  5732.4  7917.4 31479.1

Total societal costs must be at least as high as total healthcare costs because they include all healthcare costs plus productivity losses. If this is not the case, something has gone wrong.

6. Inspect the total cost distributions

We compare healthcare costs between the intervention and control groups.

Code

boxplot(
  total_costs_healthcare ~ group,
  data = df.knee,
  ylab = "Total healthcare costs (2026 CHF)",
  xlab = "Group",
  main = "Healthcare perspective"
)

The boxplot shows the centre, spread and potential outliers in each group. The intervention group has additional intervention costs, but these may be partly offset by lower use of other healthcare resources.

Finally, we examine costs from the societal perspective.

Code

boxplot(
  total_costs_societal ~ group,
  data = df.knee,
  ylab = "Total societal costs (2026 CHF)",
  xlab = "Group",
  main = "Societal perspective"
)

The comparison may change when productivity losses are included. An intervention can be more expensive from the healthcare perspective but less expensive from the societal perspective if it reduces absence or impairment during productive activities.

Costs alone do not determine cost-effectiveness

A lower total cost does not automatically make an intervention cost-effective, just as a cheaper umbrella is not necessarily a good purchase if it leaks. We still need to compare the costs with the health outcomes achieved by each group.

2.4 Clinical outcomes

Costs tell us how many resources were used, but an economic evaluation must also consider what was achieved with those resources. In this exercise, we assess outcomes in two complementary ways:

KOOS-ADL measures changes in knee-related functioning;
QALYs combine health-related quality of life with the time spent in that state of health.

KOOS-ADL ranges from 0 to 100, with higher scores indicating better functioning in activities of daily living. For each participant, clinical improvement over the complete follow-up period is calculated as:

\[ \Delta KOOS_i = KOOS_{i,24} - KOOS_{i,0}, \]

where:

$\Delta KOOS_i$ is the change in KOOS-ADL for participant $i$;
$KOOS_{i,0}$ is the baseline score;
$KOOS_{i,24}$ is the score after 24 months.

A positive value indicates improved functioning, while a negative value indicates deterioration. The difference in mean improvement between the groups will later be used to calculate the incremental cost per additional KOOS-ADL point.

The EQ-5D-5L health-state codes recorded at baseline, 6, 12 and 24 months must first be converted into utility values using an appropriate value set. A utility of 1 represents full health, 0 represents a health state equivalent to death, and negative values represent health states valued as worse than death.

QALYs are then calculated using the area-under-the-curve method. We assume that utility changes linearly between two measurement points:

\[ QALY_{t_1,t_2} = \frac{U_{t_1}+U_{t_2}}{2} \times (t_2-t_1), \]

where:

$U_{t_1}$ and $U_{t_2}$ are the utility values at two consecutive measurement points;
$(t_2-t_1)$ is the time between the measurements, expressed in years.

For example, the period from baseline to 6 months represents half a year:

\[ QALY_{0-6} = \frac{U_0+U_6}{2} \times 0.5. \]

The same calculation is used for the periods from 6 to 12 months and from 12 to 24 months.

QALYs accrued during the second year occur in the future from the perspective of the baseline treatment decision. We therefore discount the QALYs accrued between 12 and 24 months using an annual discount rate of 3%:

\[ QALY_{\text{total}} = QALY_{0-6} + QALY_{6-12} + \frac{QALY_{12-24}}{1.03}. \]

Baseline adjustment

In this exercise, we compare mean QALYs directly between the randomised groups. In a full trial-based economic evaluation, incremental QALYs would usually be adjusted for differences in baseline utility using an appropriate regression model.

The KOOS-ADL change score is treated as a clinical endpoint and is not discounted in this exercise. This is a simplification: methodological guidelines may differ in how future clinical outcomes should be handled (indeed, there is a debate about this).

Exercise

1. Calculate the change in KOOS-ADL

For each participant, calculate the change in KOOS-ADL from baseline to 24 months:

\[ \Delta KOOS_i = KOOS_{i,24} - KOOS_{i,0}. \]

Store the result in a new variable named:

koos_adl_change.

Create a boxplot comparing the intervention and control groups. You can adapt the following code:

Code

boxplot(
  koos_adl_change ~ group,
  data = df.knee,
  ylab = "Change in KOOS-ADL",
  xlab = "Group"
)

Consider whether improvement appears to differ between the two groups and whether there is variation in treatment response.

2. Convert EQ-5D-5L health states into utilities

The variables eq5d_0, eq5d_6, eq5d_12 and eq5d_24 contain five-digit EQ-5D-5L health-state codes. Use the eq5d() function from the eq5d package to convert these codes into utility values.

There is currently no Swiss EQ-5D-5L value set available in the EuroQol value-set collection. For this exercise, use the German EQ-5D-5L value set as a transparent teaching assumption.

Install and load the package:

Code

install.packages("eq5d")
library(eq5d)

The baseline utilities can be calculated as follows:

Code

df.knee$utility_0 <- eq5d(
  scores = df.knee$eq5d_0,
  version = "5L",
  type = "VT", # VT: a value set developed directly for EQ-5D-5L health states
  country = "Germany"
)

Adapt this code to calculate and store utility values for all four measurement points:

utility_0;
utility_6;
utility_12;
utility_24.

Use summary() to inspect the four utility variables. Check that the values are plausible. A utility of 1 represents full health, while values below 0 represent health states valued as worse than death.

Choice of value set

EQ-5D utilities reflect preferences from a population and are therefore influenced by the selected value set. Using the German value set is a simplification for this exercise, not an assumption that Swiss and German preferences are identical. In a real evaluation, the choice should be justified and examined in a sensitivity analysis where appropriate.

3. Calculate QALYs separately for the following periods:

baseline to 6 months;
6 to 12 months;
12 to 24 months.

Use the area-under-the-curve method:

\[ QALY_{t_1,t_2} = \frac{U_{t_1}+U_{t_2}}{2} \times (t_2-t_1). \]

Add the following variables to df.knee:

qaly_0_6;
qaly_6_12;
qaly_12_24.

For example, QALYs from baseline to 6 months can be calculated as:

Code

df.knee$qaly_0_6 <-
  ((df.knee$utility_0 + df.knee$utility_6) / 2) *
  0.5

Adapt this code for the other two periods. Pay attention to the length of each interval: the first two periods last half a year each, whereas the final period lasts one complete year.

4. Discount QALYs accrued during the second year

Use an annual discount rate of 3%. QALYs accrued during the first year remain undiscounted, while QALYs accrued between 12 and 24 months are discounted by one year.

Create:

qaly_12_24_discounted.

You can use:

Code

discount_rate <- 0.03

df.knee$qaly_12_24_discounted <-
  df.knee$qaly_12_24 /
  (1 + discount_rate)

Keep the original qaly_12_24 variable so that both the discounted and undiscounted values remain available.

5. Calculate total QALYs

Calculate the total discounted QALYs accumulated by each participant over 24 months and store the result as:

total_qalys.

The calculation is:

\[ QALY_{\text{total}} = QALY_{0-6} + QALY_{6-12} + PV(QALY_{12-24}). \]

You can start with:

Code

df.knee$total_qalys <-
  df.knee$qaly_0_6 +
  df.knee$qaly_6_12 +
  # Add the discounted QALYs from 12 to 24 months

6. Inspect the outcome distributions

Create a boxplot comparing total QALYs and KOOS-scores between the two groups. Use and adapt:

Code

boxplot(
  total_qalys ~ group,
  data = df.knee,
  ylab = "Total QALYs over 24 months",
  xlab = "Group"
)

When interpreting the results, consider:

whether the intervention group shows greater improvement in KOOS-ADL;
whether total QALYs differ between the groups;
whether participants vary substantially within each group;
whether the conclusions appear similar for the disease-specific and generic outcomes.

1. Calculate the change in KOOS-ADL

We calculate each participant’s change from baseline to 24 months. Because higher KOOS-ADL scores indicate better functioning, positive change scores represent improvement.

Code

df.knee$koos_adl_change <-
  df.knee$koos_adl_24 -
  df.knee$koos_adl_0

2. Convert EQ-5D-5L health states into utilities

We use the German EQ-5D-5L value set to convert each five-digit health-state code into a utility value.

Code

library(eq5d)

df.knee$utility_0 <- eq5d(
  scores = df.knee$eq5d_0,
  version = "5L",
  type = "VT",
  country = "Germany"
)

df.knee$utility_6 <- eq5d(
  scores = df.knee$eq5d_6,
  version = "5L",
  type = "VT",
  country = "Germany"
)

df.knee$utility_12 <- eq5d(
  scores = df.knee$eq5d_12,
  version = "5L",
  type = "VT",
  country = "Germany"
)

df.knee$utility_24 <- eq5d(
  scores = df.knee$eq5d_24,
  version = "5L",
  type = "VT",
  country = "Germany"
)

We can inspect the calculated utility values:

Code

summary(df.knee[c(
  "utility_0",
  "utility_6",
  "utility_12",
  "utility_24"
)])

   utility_0         utility_6         utility_12       utility_24     
 Min.   :-0.3730   Min.   :-0.0980   Min.   :0.0010   Min.   :-0.2930  
 1st Qu.: 0.5425   1st Qu.: 0.6790   1st Qu.:0.7448   1st Qu.: 0.7292  
 Median : 0.7625   Median : 0.8185   Median :0.8510   Median : 0.8300  
 Mean   : 0.6790   Mean   : 0.7601   Mean   :0.7975   Mean   : 0.7517  
 3rd Qu.: 0.8442   3rd Qu.: 0.9070   3rd Qu.:0.9350   3rd Qu.: 0.8940  
 Max.   : 1.0000   Max.   : 1.0000   Max.   :1.0000   Max.   : 1.0000

Code

with(df.knee, boxplot(utility_0, utility_6, utility_12, utility_24))

This output serves as a small plausibility-check: no value is > 1.

3. Calculate QALYs for each follow-up period

We calculate the area under the utility curve separately for each interval. The first two intervals last half a year each.

Code

df.knee$qaly_0_6 <-
  ((df.knee$utility_0 + df.knee$utility_6) / 2) *
  0.5

df.knee$qaly_6_12 <-
  ((df.knee$utility_6 + df.knee$utility_12) / 2) *
  0.5

The utility values at the beginning and end of each interval are averaged and multiplied by the duration of the interval. This assumes that utility changes linearly between measurement points—which is a useful approximation, although health rarely travels in a perfectly straight line.

The period from 12 to 24 months lasts one complete year:

Code

df.knee$qaly_12_24 <-
  ((df.knee$utility_12 + df.knee$utility_24) / 2) *
  1

We can inspect the QALYs accumulated during each period:

Code

summary(df.knee[c(
  "qaly_0_6",
  "qaly_6_12",
  "qaly_12_24"
)])

    qaly_0_6         qaly_6_12         qaly_12_24     
 Min.   :-0.1177   Min.   :0.06725   Min.   :-0.1295  
 1st Qu.: 0.3161   1st Qu.:0.34187   1st Qu.: 0.7208  
 Median : 0.3850   Median :0.41413   Median : 0.8267  
 Mean   : 0.3598   Mean   :0.38942   Mean   : 0.7746  
 3rd Qu.: 0.4290   3rd Qu.:0.45025   3rd Qu.: 0.9075  
 Max.   : 0.5000   Max.   :0.50000   Max.   : 1.0000

4. Discount QALYs accrued during the second year

We apply an annual discount rate of 3% to the QALYs accrued between 12 and 24 months.

Code

discount_rate <- 0.03

df.knee$qaly_12_24_discounted <-
  df.knee$qaly_12_24 /
  (1 + discount_rate)

The discounted values must be slightly lower than the undiscounted year-2 QALYs. The effect is modest because only one year is discounted, but over longer time horizons the difference becomes increasingly important.

5. Calculate total QALYs

Total QALYs consist of the undiscounted QALYs accrued during the first year and the discounted QALYs accrued during the second year.

Code

df.knee$total_qalys <-
  df.knee$qaly_0_6 +
  df.knee$qaly_6_12 +
  df.knee$qaly_12_24_discounted

A participant who remained in full health throughout the complete follow-up would accumulate slightly less than two discounted QALYs:

\[ 1 + \frac{1}{1.03} \approx 1.971. \]

Code

max(df.knee$total_qalys)

[1] 1.970874

The observed values should generally be lower because the participants experience some limitations in health-related quality of life.

6. Inspect the outcome distributions

Code

boxplot(
  koos_adl_change ~ group,
  data = df.knee,
  ylab = "Change in KOOS-ADL",
  xlab = "Group",
  main = "Change in knee-related functioning"
)

The intervention group tends to show greater improvement in KOOS-ADL. However, the overlap between the groups indicates that participants differ in their response to treatment.

Code

boxplot(
  total_qalys ~ group,
  data = df.knee,
  ylab = "Total QALYs over 24 months",
  xlab = "Group",
  main = "QALYs by treatment group"
)

The intervention group tends to accumulate more QALYs, although the distributions overlap.

Clinical improvement and QALY gain are related but not identical

A participant may improve in KOOS-ADL without showing an equally large gain in QALYs. EQ-5D-5L measures broader aspects of health and may be less sensitive to specific changes in knee-related functioning.

This is why we retain both outcomes: one keeps a close eye on the knee, while the other asks how life as a whole is going.

Still, there is quite strong correlation between the two:

Code

cor.test(df.knee$koos_adl_24, df.knee$utility_24)


    Pearson's product-moment correlation

data:  df.knee$koos_adl_24 and df.knee$utility_24
t = 15.729, df = 198, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.6765089 0.8011919
sample estimates:
      cor 
0.7452967

2.5 Incremental cost-effectiveness

Full health economic evaluation compares the additional costs and outcomes of the intervention with those of the control strategy. We therefore focus on differences between the group means.

Incremental costs are calculated as:

\[ \Delta C = \bar{C}_{\text{intervention}} - \bar{C}_{\text{control}}, \]

and incremental effects as:

\[ \Delta E = \bar{E}_{\text{intervention}} - \bar{E}_{\text{control}}, \]

where:

$\Delta C$ is the difference in mean costs;
$\Delta E$ is the difference in mean effects;
$\bar{C}$ and $\bar{E}$ are the corresponding group means.

Using KOOS-ADL as the outcome, the incremental cost-effectiveness ratio is:

\[ ICER = \frac{\Delta C}{\Delta KOOS}. \]

The ICER represents the additional cost per additional KOOS-ADL point achieved by the intervention.

Using QALYs as the outcome, the incremental cost-utility ratio is:

\[ ICUR = \frac{\Delta C}{\Delta QALY}. \]

The ICUR represents the additional cost per additional QALY gained. Both ratios can be calculated from the healthcare and societal perspectives.

Interpret the differences before the ratio

A negative ICER or ICUR is not self-explanatory. It can indicate that the intervention is less costly and more effective, but it can also indicate that it is more costly and less effective.

Always inspect the signs of the incremental costs and effects before interpreting the ratio. Ratios are good at division but surprisingly poor at explaining themselves.

Bootstrapping

The ICER and ICUR are point estimates based on the observed trial sample. To explore how much these estimates might vary across different samples, participants can be repeatedly resampled with replacement using bootstrapping.

Non-parametric bootstrapping repeatedly draws samples of participants with replacement from the observed data, without assuming a particular probability distribution for costs or effects. Each selected participant contributes their complete set of costs and outcomes, and the incremental costs and effects are recalculated in every bootstrap sample. Repeating this process many times approximates the sampling distribution of the economic results.

For each treatment group $g$, participants are selected as:

\[ I_{gk}^{(b)} \sim \operatorname{Uniform}\{1,\ldots,n_g\}, \qquad X_{gk}^{*(b)} = X_{g,I_{gk}^{(b)}}, \]

where:

$I_{gk}^{(b)}$ is the participant index selected from group $g$ on draw $k$ in bootstrap sample $b$;
$n_g$ is the number of participants in group $g$;
$X_{gk}^{*(b)}$ is the complete record of the selected participant.

Each participant within a group has the same probability of being selected on every draw:

\[ \Pr\left(I_{gk}^{(b)}=i\right) = \frac{1}{n_g}. \]

In this trial, each group contains 100 participants. The probability of selecting a particular participant on one draw within their group is therefore:

\[ \frac{1}{100} = 0.01. \]

For every bootstrap sample, we make 100 draws from the intervention group and 100 draws from the control group. Sampling is performed with replacement, so the same participant may be selected several times, while another participant may not be selected at all.

Why resample within each group?

Resampling separately within the intervention and control groups preserves the original group sizes. Every selected participant retains their complete set of costs and outcomes, maintaining the observed relationships between these variables.

Cost-effectiveness plane

The resulting estimates of each bootstrap sample can be displayed on a cost-effectiveness plane:

the horizontal axis shows incremental effects;
the vertical axis shows incremental costs.

Quadrant	Incremental costs	Incremental effects	Interpretation
North-east	Higher	Better	More effective and more costly
South-east	Lower	Better	Dominant
South-west	Lower	Worse	Less effective and less costly
North-west	Higher	Worse	Dominated

The location and spread of the bootstrap estimates show both the likely economic result and the uncertainty surrounding it.

When an intervention is more effective but also more costly, its cost-effectiveness depends on how much a decision-maker is willing to pay for the additional benefit.

For QALYs, the intervention is considered cost-effective when:

\[ \lambda \times \Delta QALY - \Delta C > 0, \]

where:

$\lambda$ is the willingness to pay for one additional QALY;
$\Delta QALY$ is the incremental QALY gain;
$\Delta C$ is the incremental cost.

The cost-effectiveness acceptability curve shows the proportion of bootstrap samples in which this condition is satisfied across a range of willingness-to-pay thresholds.

The curve therefore shows the probability that the intervention is cost-effective at each threshold. It does not prove that the intervention is cost-effective; it shows how confident we can be under different assumptions about the value of an additional QALY.

Exercise

In this section, you will calculate the incremental costs and effects, derive the ICER and ICUR, and use bootstrapping to examine how the economic results might vary across repeated samples.

1. Calculate the mean costs and outcomes in each group

Calculate the mean values for the intervention and control groups separately for:

total_costs_healthcare;
total_costs_societal;
koos_adl_change;
total_qalys.

Although cost data are often right-skewed, use the mean rather than the median (!). Economic evaluation is concerned with the expected cost per participant, and total costs must balance across the population—even if a few hospital admissions try their best to dominate the calculation :-)

For example, mean healthcare costs in the intervention group can be calculated as:

Code

mean_cost_healthcare_intervention <- mean(
  df.knee$total_costs_healthcare[
    df.knee$group == "Intervention"
  ]
)

Adapt this code to calculate the corresponding mean for the control group and for the other costs and outcomes.

2. Calculate incremental costs and effects

Calculate the intervention-group mean minus the control-group mean for:

healthcare costs;
societal costs;
KOOS-ADL change;
QALYs.

Store the results as:

delta_cost_healthcare;
delta_cost_societal;
delta_koos;
delta_qaly.

For example:

Code

delta_cost_healthcare <-
  mean_cost_healthcare_intervention -
  mean_cost_healthcare_control

Interpret the sign of every result:

a positive incremental cost means that the intervention is more costly;
a positive incremental effect means that it is more effective;

3. Calculate the ICER and ICUR

Calculate the ICER using KOOS-ADL change as the outcome:

\[ ICER = \frac{\Delta C}{\Delta KOOS}. \]

Calculate one ICER from each perspective and store the results as:

icer_healthcare;
icer_societal.

For example:

Code

icer_healthcare <-
  delta_cost_healthcare /
  delta_koos

Next, calculate the ICUR using QALYs as the outcome:

\[ ICUR = \frac{\Delta C}{\Delta QALY}. \]

Store the results as:

icur_healthcare;
icur_societal.

Before interpreting any ratio, examine the signs of its numerator and denominator. State whether the intervention is:

more effective and more costly;
more effective and less costly;
less effective and less costly;
less effective and more costly.

If an incremental effect is close to zero, the corresponding ratio may become very large and unstable. That is mathematics being technically correct but not especially helpful.

4. Prepare the data for bootstrapping

Create separate datasets for the intervention and control groups:

Code

df.intervention <- df.knee[
  df.knee$group == "Intervention",
]

df.control <- df.knee[
  df.knee$group == "Control",
]

Choose the number of bootstrap samples and create an empty dataset in which the results can be stored:

Code

set.seed(2026) # makes the bootstrap reproducible

n_boot <- 2000

bootstrap.results <- data.frame(
  delta_cost_healthcare = numeric(n_boot),
  delta_cost_societal = numeric(n_boot),
  delta_koos = numeric(n_boot),
  delta_qaly = numeric(n_boot)
)

Using set.seed() ensures that the same bootstrap samples can be reproduced. We use 2000 repetitions for this exercise. More repetitions may produce more stable estimates but also require more computing time—and more patience from the laptop (not so much for the new mac generation…).

5. Generate the bootstrap samples

The following code performs the complete bootstrap procedure for the healthcare perspective. In every repetition, it:

samples participants with replacement within each treatment group;
creates new intervention and control samples;
calculates incremental healthcare costs;
calculates incremental KOOS-ADL and QALYs;
stores the results.

Code

set.seed(2026)

n_boot <- 2000

bootstrap.results <- data.frame(
  delta_cost_healthcare = numeric(n_boot),
  delta_cost_societal = numeric(n_boot),
  delta_koos = numeric(n_boot),
  delta_qaly = numeric(n_boot)
)

for (b in 1:n_boot) {
  
  index_intervention <- sample(
    1:nrow(df.intervention),
    size = nrow(df.intervention),
    replace = TRUE
  )
  
  index_control <- sample(
    1:nrow(df.control),
    size = nrow(df.control),
    replace = TRUE
  )
  
  bootstrap.intervention <-
    df.intervention[index_intervention, ]
  
  bootstrap.control <-
    df.control[index_control, ]
  
  bootstrap.results$delta_cost_healthcare[b] <-
    mean(bootstrap.intervention$total_costs_healthcare) -
    mean(bootstrap.control$total_costs_healthcare)
  
  bootstrap.results$delta_koos[b] <-
    mean(bootstrap.intervention$koos_adl_change) -
    mean(bootstrap.control$koos_adl_change)
  
  bootstrap.results$delta_qaly[b] <-
    mean(bootstrap.intervention$total_qalys) -
    mean(bootstrap.control$total_qalys)
}

Adapt the code by adding one calculation inside the for loop for:

bootstrap.results$delta_cost_societal[b].

Use total_costs_societal instead of total_costs_healthcare. The same resampled participants must be used for both perspectives.

Add the calculation after the healthcare-cost calculation:

Code

bootstrap.results$delta_cost_societal[b] <-
  # Mean societal costs in the bootstrap intervention group
  # minus mean societal costs in the bootstrap control group

Only the cost variable changes. The incremental KOOS-ADL and QALY estimates are the same for both perspectives because the analytical perspective determines which costs are included, not how health outcomes are measured.

Use the same bootstrap samples

Do not run a separate bootstrap for the societal perspective. Healthcare and societal costs must be calculated using the same resampled participants. Otherwise, differences between the perspectives would partly reflect different random samples rather than the inclusion of productivity costs.

6. Inspect the bootstrap results

Use summary() to inspect the bootstrap estimates:

Code

summary(bootstrap.results)

Compare the centre of each bootstrap distribution with the corresponding point estimate from the original sample. They should be reasonably similar, although they will not be exactly identical.

You can also calculate percentile-based 95% bootstrap intervals. For example:

Code

quantile(
  bootstrap.results$delta_cost_healthcare,
  probs = c(0.025, 0.975)
)

Adapt this code to inspect the intervals for incremental societal costs, KOOS-ADL and QALYs.

7. Create a cost-effectiveness plane

Create a cost-effectiveness plane using incremental QALYs on the horizontal axis and incremental healthcare costs on the vertical axis:

Code

plot(
  bootstrap.results$delta_qaly,
  bootstrap.results$delta_cost_healthcare,
  pch = 16,
  col = rgb(0, 0, 1, 0.20),
  xlab = "Incremental QALYs",
  ylab = "Incremental healthcare costs (2026 CHF)"
)

abline(h = 0, v = 0, lty = 2)

The transparent colour helps reveal where the bootstrap estimates are concentrated. Add the point estimate from the original sample:

Code

points(
  delta_qaly,
  delta_cost_healthcare,
  pch = 19,
  col = "red"
)

Adapt the plot to create a second cost-effectiveness plane using delta_cost_societal.

Interpret:

the quadrant containing most bootstrap estimates;
the amount of spread in incremental costs and QALYs;
whether estimates occur in several quadrants;
whether the conclusion appears more favourable from one perspective than the other.

Of course, you could do the same for the KOOS-ADL.

8. Calculate the probability of cost-effectiveness

Create a sequence of willingness-to-pay thresholds ranging from CHF 0 to CHF 200,000 per QALY:

Code

wtp <- seq(
  from = 0,
  to = 200000,
  by = 5000
)

At each threshold, the intervention is considered cost-effective in a bootstrap sample when:

\[ \lambda \times \Delta QALY - \Delta C > 0. \]

The following code calculates the probability of cost-effectiveness from the healthcare perspective:

Code

prob_ce_healthcare <- numeric(length(wtp))

for (j in 1:length(wtp)) {
  
  prob_ce_healthcare[j] <- mean(
    wtp[j] * bootstrap.results$delta_qaly -
      bootstrap.results$delta_cost_healthcare > 0
  )
}

The expression inside mean() returns TRUE when the intervention is cost-effective and FALSE when it is not. In R, TRUE is treated as 1 and FALSE as 0. The mean therefore gives the proportion of bootstrap samples in which the intervention is cost-effective.

For example, a result of 0.72 means that the intervention was cost-effective in 72% of the bootstrap samples at that willingness-to-pay threshold.

Now adapt the code to calculate the probability from the societal perspective:

create an empty vector named prob_ce_societal;
use another for loop over the same willingness-to-pay thresholds;
replace delta_cost_healthcare with delta_cost_societal;
keep delta_qaly unchanged.

You can start with:

Code

prob_ce_societal <- numeric(length(wtp))

for (j in 1:length(wtp)) {
  
  prob_ce_societal[j] <- mean(
    # Add the cost-effectiveness condition
    # using incremental societal costs
  )
}

Only the cost variable changes between the two perspectives. The willingness-to-pay thresholds and incremental QALYs remain the same. Thankfully, the bootstrap has already done the hard work; we are now merely asking it the same question from a different perspective.

9. Create the cost-effectiveness acceptability curve

Plot the probability of cost-effectiveness against the willingness-to-pay threshold:

Code

plot(
  wtp,
  prob_ce_healthcare,
  type = "l",
  lwd = 2,
  ylim = c(0, 1),
  xlab = "Willingness to pay per QALY (CHF)",
  ylab = "Probability cost-effective"
)

Do the same for the societal perspective. Use the two graphs to answer the following questions:

What is the probability that the intervention is cost-effective when the willingness to pay is CHF 0 per QALY?
How does this probability change as willingness to pay increases?
Does the conclusion differ between the healthcare and societal perspectives?
Does the curve reach a probability close to 1, or does substantial uncertainty remain?
At which willingness-to-pay threshold does the probability first exceed 50%, if it does?

Remember that the CEAC describes uncertainty rather than making the decision itself. A probability of 70% does not mean that 70% of patients benefit; it means that 70% of the bootstrap samples classify the intervention as cost-effective at that threshold.

1. Calculate the mean costs and outcomes in each group

We first calculate the mean costs and outcomes separately for the intervention and control groups.

Code

mean_cost_healthcare_intervention <- mean(
  df.knee$total_costs_healthcare[
    df.knee$group == "Intervention"
  ]
)

mean_cost_healthcare_control <- mean(
  df.knee$total_costs_healthcare[
    df.knee$group == "Control"
  ]
)

mean_cost_societal_intervention <- mean(
  df.knee$total_costs_societal[
    df.knee$group == "Intervention"
  ]
)

mean_cost_societal_control <- mean(
  df.knee$total_costs_societal[
    df.knee$group == "Control"
  ]
)

These values represent the expected healthcare and societal costs per participant in each group.

Next, we calculate the mean KOOS-ADL change:

Code

mean_koos_intervention <- mean(
  df.knee$koos_adl_change[
    df.knee$group == "Intervention"
  ]
)

mean_koos_control <- mean(
  df.knee$koos_adl_change[
    df.knee$group == "Control"
  ]
)

A higher mean change indicates a greater average improvement in knee-related functioning.

Finally, we calculate the mean QALYs:

Code

mean_qaly_intervention <- mean(
  df.knee$total_qalys[
    df.knee$group == "Intervention"
  ]
)

mean_qaly_control <- mean(
  df.knee$total_qalys[
    df.knee$group == "Control"
  ]
)

Higher mean QALYs indicate that participants experienced better health-related quality of life over the 24-month follow-up period.

We can present all group means together (not needed for further calculations):

Code

group.means <- data.frame(
  group = c("Intervention", "Control"),
  healthcare_costs = c(
    mean_cost_healthcare_intervention,
    mean_cost_healthcare_control
  ),
  societal_costs = c(
    mean_cost_societal_intervention,
    mean_cost_societal_control
  ),
  koos_change = c(
    mean_koos_intervention,
    mean_koos_control
  ),
  qalys = c(
    mean_qaly_intervention,
    mean_qaly_control
  )
)

group.means

         group healthcare_costs societal_costs koos_change    qalys
1 Intervention         2337.469       4489.769      10.539 1.608009
2      Control         2771.742       6975.014       3.166 1.394478

The table allows us to compare the costs and outcomes directly. However, economic evaluation is interested in the differences between the groups, which we calculate next.

2. Calculate incremental costs and effects

Incremental values are calculated as the intervention-group mean minus the control-group mean.

Code

delta_cost_healthcare <-
  mean_cost_healthcare_intervention -
  mean_cost_healthcare_control

delta_cost_societal <-
  mean_cost_societal_intervention -
  mean_cost_societal_control

delta_koos <-
  mean_koos_intervention -
  mean_koos_control

delta_qaly <-
  mean_qaly_intervention -
  mean_qaly_control

We inspect the four incremental estimates together:

Code

incremental.results <- data.frame(
  outcome = c(
    "Healthcare costs",
    "Societal costs",
    "KOOS-ADL change",
    "QALYs"
  ),
  difference = c(
    delta_cost_healthcare,
    delta_cost_societal,
    delta_koos,
    delta_qaly
  )
)

incremental.results

           outcome    difference
1 Healthcare costs  -434.2730097
2   Societal costs -2485.2446047
3  KOOS-ADL change     7.3730000
4            QALYs     0.2135307

Positive incremental effects indicate that the intervention produces better outcomes on average. Positive incremental costs indicate that it is more expensive, whereas negative incremental costs indicate cost savings. In this case, all four estimates are in favour of the intervention group. However, note that the same intervention may appear more or less costly depending on the perspective because productivity losses are included only in societal costs. Perspective matters: costs do not disappear simply because the healthcare budget cannot see them.

3. Calculate the ICER and ICUR

The ICER expresses the incremental cost per additional KOOS-ADL point.

Code

icer_healthcare <-
  delta_cost_healthcare /
  delta_koos

icer_societal <-
  delta_cost_societal /
  delta_koos

icer.results <- data.frame(
  perspective = c("Healthcare", "Societal"),
  icer = c(
    icer_healthcare,
    icer_societal
  )
)

icer.results

  perspective       icer
1  Healthcare  -58.90045
2    Societal -337.07373

The healthcare ICER uses only healthcare costs, whereas the societal ICER additionally considers productivity losses. The denominator is identical because changing perspective affects the included costs, not the measured clinical effect.

The ICUR expresses the incremental cost per additional QALY.

Code

icur_healthcare <-
  delta_cost_healthcare /
  delta_qaly

icur_societal <-
  delta_cost_societal /
  delta_qaly

icur.results <- data.frame(
  perspective = c("Healthcare", "Societal"),
  icur = c(
    icur_healthcare,
    icur_societal
  )
)

icur.results

  perspective       icur
1  Healthcare  -2033.773
2    Societal -11638.817

Before interpreting these ratios, we inspect the signs of the incremental costs and effects:

positive costs and positive effects: the intervention is more effective and more costly;
negative costs and positive effects: the intervention is dominant (which is the case here);
negative costs and negative effects: the intervention is less effective and less costly;
positive costs and negative effects: the intervention is dominated.

A negative ratio must therefore never be interpreted without examining its numerator and denominator. The minus sign tells us that the directions differ, but it does not tell us whether that is good news or bad news.

4. Prepare the data for bootstrapping

Ok - let’s get the nerdy part done. We create separate datasets for the two treatment groups. Bootstrapping will then be performed separately within each group.

Code

df.intervention <- df.knee[
  df.knee$group == "Intervention",
]

df.control <- df.knee[
  df.knee$group == "Control",
]

nrow(df.intervention)

[1] 100

Code

nrow(df.control)

[1] 100

Both groups contain 100 participants. Resampling separately preserves these original group sizes.

We specify the number of bootstrap repetitions and prepare a dataset for storing the results:

Code

set.seed(2026)

n_boot <- 2000

bootstrap.results <- data.frame(
  delta_cost_healthcare = numeric(n_boot),
  delta_cost_societal = numeric(n_boot),
  delta_koos = numeric(n_boot),
  delta_qaly = numeric(n_boot)
)

At this point, all entries are zero. They will be replaced with the results from the 2,000 bootstrap samples.

5. Generate the bootstrap samples

In each repetition, we resample participants with replacement within their treatment group. Each selected participant retains their complete set of costs and outcomes.

Code

set.seed(2026)

for (b in 1:n_boot) {
  
  index_intervention <- sample(
    1:nrow(df.intervention),
    size = nrow(df.intervention),
    replace = TRUE
  )
  
  index_control <- sample(
    1:nrow(df.control),
    size = nrow(df.control),
    replace = TRUE
  )
  
  bootstrap.intervention <-
    df.intervention[index_intervention, ]
  
  bootstrap.control <-
    df.control[index_control, ]
  
  bootstrap.results$delta_cost_healthcare[b] <-
    mean(bootstrap.intervention$total_costs_healthcare) -
    mean(bootstrap.control$total_costs_healthcare)
  
  bootstrap.results$delta_cost_societal[b] <- # This is the new part
    mean(bootstrap.intervention$total_costs_societal) -
    mean(bootstrap.control$total_costs_societal)
  
  bootstrap.results$delta_koos[b] <-
    mean(bootstrap.intervention$koos_adl_change) -
    mean(bootstrap.control$koos_adl_change)
  
  bootstrap.results$delta_qaly[b] <-
    mean(bootstrap.intervention$total_qalys) -
    mean(bootstrap.control$total_qalys)
}

The same resampled participants are used for healthcare costs, societal costs, KOOS-ADL and QALYs. This preserves the relationships between costs and outcomes within participants.

6. Inspect the bootstrap results

We first inspect the distributions of the bootstrap estimates:

Code

summary(bootstrap.results)

 delta_cost_healthcare delta_cost_societal   delta_koos       delta_qaly     
 Min.   :-2063.6       Min.   :-5106.3     Min.   : 4.626   Min.   :0.04305  
 1st Qu.: -701.9       1st Qu.:-2966.3     1st Qu.: 6.795   1st Qu.:0.18305  
 Median : -432.4       Median :-2477.4     Median : 7.362   Median :0.21334  
 Mean   : -436.4       Mean   :-2486.6     Mean   : 7.371   Mean   :0.21409  
 3rd Qu.: -173.6       3rd Qu.:-1996.6     3rd Qu.: 7.958   3rd Qu.:0.24425  
 Max.   :  871.0       Max.   : -258.5     Max.   :10.178   Max.   :0.38472

The means and medians of the bootstrap distributions should be reasonably close to the point estimates calculated from the original trial sample (which is the case, YEAAHH).

Code

delta_cost_healthcare

[1] -434.273

Code

delta_cost_societal

[1] -2485.245

Code

delta_koos

[1] 7.373

Code

delta_qaly

[1] 0.2135307

Small differences are expected because the bootstrap samples are generated randomly.

We calculate percentile-based 95% bootstrap intervals using the 2.5th and 97.5th percentiles:

Code

quantile(
  bootstrap.results$delta_cost_healthcare,
  probs = c(0.025, 0.975)
)

      2.5%      97.5% 
-1213.5465   345.9892

Code

quantile(
  bootstrap.results$delta_cost_societal,
  probs = c(0.025, 0.975)
)

     2.5%     97.5% 
-3902.664 -1112.079

Code

quantile(
  bootstrap.results$delta_koos,
  probs = c(0.025, 0.975)
)

    2.5%    97.5% 
5.631975 9.050500

Code

quantile(
  bootstrap.results$delta_qaly,
  probs = c(0.025, 0.975)
)

     2.5%     97.5% 
0.1266750 0.3049837

Wide intervals indicate substantial sampling uncertainty. If an interval includes zero, the bootstrap results include both positive and negative differences. This does not automatically mean that there is “no effect”; it means that the direction of the estimate is uncertain given the available sample. We see that for health care costs, the interval includes 0.

7. Create cost-effectiveness planes

We first create the cost-effectiveness plane from the healthcare perspective.

Code

plot(
  bootstrap.results$delta_qaly,
  bootstrap.results$delta_cost_healthcare,
  pch = 16,
  col = rgb(0, 0, 1, 0.20),
  xlab = "Incremental QALYs",
  ylab = "Incremental healthcare costs (2026 CHF)",
  main = "Healthcare perspective",
  ylim = c(-2000, 2000),
  xlim = c(-0.05, 0.4)
)

abline(
  h = 0,
  v = 0,
  lty = 2
)

points(
  delta_qaly,
  delta_cost_healthcare,
  pch = 19,
  col = "red"
)

Each blue point represents the incremental costs and QALYs from one bootstrap sample. The red point represents the result from the original trial sample.

We then create the corresponding plane from the societal perspective:

Code

plot(
  bootstrap.results$delta_qaly,
  bootstrap.results$delta_cost_societal,
  pch = 16,
  col = rgb(0, 0.5, 0, 0.20),
  xlab = "Incremental QALYs",
  ylab = "Incremental societal costs (2026 CHF)",
  main = "Societal perspective",
  xlim = c(-0.05, 0.4),
  ylim = c(-5000, 5000)
)

abline(
  h = 0,
  v = 0,
  lty = 2
)

points(
  delta_qaly,
  delta_cost_societal,
  pch = 19,
  col = "red"
)

The horizontal position of the bootstrap estimates is identical in both figures because both use incremental QALYs. Their vertical positions differ because the healthcare and societal perspectives include different costs.

The plots should be interpreted according to the quadrants:

north-east: more effective and more costly;
south-east: more effective and less costly (our case);
south-west: less effective and less costly;
north-west: less effective and more costly.

A concentrated cloud indicates relatively little sampling uncertainty. A cloud spread across several quadrants indicates that the economic conclusion is less certain. The bootstrap cloud is therefore considerably more informative than one point estimate standing alone.

8. Calculate the probability of cost-effectiveness

We define willingness-to-pay thresholds between CHF 0 and CHF 200,000 per QALY.

Code

wtp <- seq(
  from = 0,
  to = 200000,
  by = 5000
)

We calculate the probability of cost-effectiveness from the healthcare perspective:

Code

prob_ce_healthcare <- numeric(length(wtp))

for (j in 1:length(wtp)) {
  
  prob_ce_healthcare[j] <- mean(
    wtp[j] * bootstrap.results$delta_qaly -
      bootstrap.results$delta_cost_healthcare > 0
  )
}

At each threshold, the condition is evaluated separately in all 2,000 bootstrap samples. Because TRUE is treated as 1 and FALSE as 0, the mean is the proportion of samples in which the intervention is cost-effective.

We repeat the calculation using societal costs:

Code

prob_ce_societal <- numeric(length(wtp))

for (j in 1:length(wtp)) {
  
  prob_ce_societal[j] <- mean(
    wtp[j] * bootstrap.results$delta_qaly -
      bootstrap.results$delta_cost_societal > 0
  )
}

Only the incremental cost variable changes. The incremental QALY estimates and willingness-to-pay thresholds remain the same.

We can inspect the calculated probabilities:

Code

ceac.results <- data.frame(
  wtp = wtp,
  probability_healthcare = prob_ce_healthcare,
  probability_societal = prob_ce_societal
)

summary(ceac.results)

      wtp         probability_healthcare probability_societal
 Min.   :     0   Min.   :0.8600         Min.   :1           
 1st Qu.: 50000   1st Qu.:1.0000         1st Qu.:1           
 Median :100000   Median :1.0000         Median :1           
 Mean   :100000   Mean   :0.9966         Mean   :1           
 3rd Qu.:150000   3rd Qu.:1.0000         3rd Qu.:1           
 Max.   :200000   Max.   :1.0000         Max.   :1

At a willingness to pay of CHF 0, the intervention is considered cost-effective only when it is less costly than the control strategy. As the threshold increases, additional QALYs receive a higher monetary value.

9. Create the cost-effectiveness acceptability curves

We first plot the CEAC from the healthcare perspective:

Code

plot(
  wtp,
  prob_ce_healthcare,
  type = "l",
  lwd = 2,
  ylim = c(0, 1),
  xlab = "Willingness to pay per QALY (CHF)",
  ylab = "Probability cost-effective",
  main = "CEAC: healthcare perspective"
)

abline(
  h = 0.5,
  lty = 2
)

The dashed horizontal line indicates a probability of 50%. The point at which the curve crosses this line shows the approximate threshold at which the intervention becomes more likely than not to be cost-effective.

We then create the CEAC from the societal perspective:

Code

plot(
  wtp,
  prob_ce_societal,
  type = "l",
  lwd = 2,
  col = "darkgreen",
  ylim = c(0, 1),
  xlab = "Willingness to pay per QALY (CHF)",
  ylab = "Probability cost-effective",
  main = "CEAC: societal perspective"
)

abline(
  h = 0.5,
  lty = 2
)

Differences between the curves arise because the societal perspective includes productivity losses. If the intervention reduces these losses, it may have a higher probability of being cost-effective from the societal perspective.

The CEAC does not show the proportion of participants who benefit from treatment. It shows the proportion of bootstrap samples in which the intervention provides positive economic value at each willingness-to-pay threshold.

The probability of cost-effectiveness is high from both perspectives across most of the willingness-to-pay range.

Point estimate versus decision uncertainty

The ICUR provides one estimate of the additional cost per QALY. The CEAC adds another layer by showing how certain we are that the intervention is cost-effective across different willingness-to-pay thresholds.

The point estimate tells us where we landed; the CEAC tells us how much the landing spot might move if the trial were repeated.

2.6 Sensitivity analysis

The results of an economic evaluation depend partly on assumptions and input values that may be uncertain. Sensitivity analysis examines whether the conclusions change when these assumptions or values are varied.

This differs from bootstrapping. Bootstrapping examines sampling uncertainty by resampling participants, whereas sensitivity analysis examines parameter and assumption uncertainty, such as uncertainty about unit costs, discount rates or the monetary value of productive time.

To keep this exercise manageable, there is no additional task on sensitivity analysis. The following sections nevertheless describe what different forms of sensitivity analysis could look like in the context of our evaluation.

One-way sensitivity analysis

In a one-way sensitivity analysis, one parameter is changed while all other parameters remain fixed at their base-case values. This approach is useful when the precise value of a particular parameter is uncertain.

For example, the base-case cost of one intervention session is CHF 55. We could repeat the analysis using:

CHF 45 per session;
CHF 65 per session;
CHF 75 per session.

The same approach could be used for other uncertain parameters, such as the discount rate, cost of a hospital admission or value of one productive hour.

After changing the selected parameter, we would recalculate the affected cost components, total costs and the ICUR. Comparing the results across the different parameter values shows how strongly the economic conclusion depends on that particular assumption.

If the conclusion remains similar across all values, the result is relatively robust to uncertainty in that parameter. If the conclusion changes substantially, the parameter is influential and should be examined particularly carefully.

I you have a proper R-script, this is quit easy: just change the paramter and run the whole script again 🤓

Multi-way sensitivity analysis

In a multi-way sensitivity analysis, several parameters are changed simultaneously. This is often presented through coherent scenarios.

For this evaluation, we could define:

a “lower-cost scenario” with lower intervention, hospital and productivity unit costs;
a “higher-cost scenario” with higher unit costs and an alternative discount rate.

The complete analysis would be repeated for each scenario. This approach shows how the results change when several plausible assumptions occur together.

The scenarios should be clinically and economically meaningful. Simply changing every value in the direction that produces the most favourable result may create a mathematically possible but unrealistic world, rather like assuming cheaper treatment, perfect adherence etc.

Probabilistic sensitivity analysis

In a probabilistic sensitivity analysis, uncertain parameters are represented by probability distributions rather than by single values. Values (i.e. the parameters) are repeatedly drawn from these distributions, and the economic evaluation is recalculated after every draw.

For example, we could assign distributions to:

the cost of an intervention session;
the cost of a hospital admission;
the monthly medication cost;
the value of one productive hour.

In each simulation, a new value for the parameter would be drawn for every uncertain parameter. The resulting distribution of incremental costs and QALYs would show the combined effect of parameter uncertainty.

Probabilistic sensitivity analysis is particularly important in model-based evaluations, such as decision trees and Markov models. In this trial-based evaluation, sampling uncertainty has already been examined using bootstrapping. Combining bootstrapping with probabilistic parameter variation would be possible, but it would add some complexity.

Three related questions

One-way sensitivity analysis: What happens if one uncertain parameter changes?
Multi-way sensitivity analysis: What happens if several assumptions change together?
Probabilistic sensitivity analysis: What range of results arises when all uncertain parameters vary simultaneously according to specified probability distributions?

Sensitivity analysis does not identify the “correct” set of assumptions. It shows whether the economic conclusion depends heavily on choices that could reasonably have been different.

3 Markov model

Trial-based economic evaluations are limited to the participants, treatments and follow-up period observed in the trial. Decision-makers may nevertheless be interested in consequences that occur over a much longer period. A treatment may, for example, delay surgery, reduce later complications or affect costs and quality of life many years after the trial has ended.

A Markov model represents the development of a population over time using a set of mutually exclusive health states. At any point in time, each person must be in exactly one state. During each model cycle, people may remain in their current state or move to another state according to predefined transition probabilities.

In a Markov model, we do not follow individual patients. Instead, we follow the proportion of a hypothetical cohort occupying each health state. Each health state can be assigned a cost and a utility. Costs and QALYs are accumulated over repeated cycles and can therefore be projected beyond the follow-up period of a clinical trial.

Markov models are particularly useful when:

a condition develops over a long period;
patients can move between clinically meaningful health states;
important events may occur repeatedly;
the timing of events affects costs or outcomes;
the consequences extend beyond the available trial data.

The Markov assumption

A conventional Markov model assumes that the probability of moving to another state depends on the current state, but not on the complete history of how the person arrived there. This is known as the “memoryless property”.

The assumption simplifies the model considerably, although patients themselves are rarely quite so willing to forget their medical history.

In this exercise, we use a simplified version of the model developed by Vetsch et al. (2023). The original evaluation examined whether implementing guideline-recommended non-surgical treatments for knee osteoarthritis would be cost-effective from the perspective of Swiss statutory healthcare.

We compare two strategies:

the current model of care, in which non-surgical treatment primarily consists of pain medication and advice;
an optimised model of care, in which structured exercise, education and other guideline-recommended non-surgical treatments are implemented.

The optimised treatment is initially more expensive but produces a larger improvement in health-related quality of life. We additionally assume that the optimised treatment delays total knee replacement by two years. In both strategies, some patients eventually undergo total knee replacement, and a small proportion subsequently require revision surgery.

The original publication used an age-specific lifetime model with a 70-year time horizon and additional tunnel states to represent delays in total knee replacement. Reproducing the complete model would require considerably more time and possibly a modest supply of emergency coffee.

For this teaching exercise, we therefore use the following model structure and assumptions:

the cohort enters the model at age 50;
one model cycle represents one year;
the model runs for 40 cycles, ending when the cohort reaches age 90;
age-specific probabilities are used for total knee replacement, revision surgery and mortality;
death is included as an absorbing state;
current care is compared with optimised care that is assumed to delay total knee replacement by two years;
costs and utilities are assigned to the health states using values from the published evaluation;
costs and QALYs are discounted at 3% annually;
probabilistic sensitivity analysis is not included.

We define the starting population as a cohort of 400 people aged 50 years who have knee osteoarthritis and are beginning non-surgical treatment.

The objective is not to reproduce the published results exactly. Instead, the exercise uses the same clinical decision and model parameters to demonstrate how a Markov model is constructed, calculated and interpreted.

We will construct the model step by step. Closely related elements are combined to avoid unnecessary repetition:

Define health states and model structure and specify transition matrices.
Follow the cohort over repeated cycles.
Assign costs and utilities to the health states.
Calculate total discounted costs and QALYs.
Compare current and optimised care.

Each step contains:

a short explanation of the relevant theory;
a task with almost-complete example code;
a solution with interpretation.

3.1 1. Define health states and model structure and specify transition matrices.

A Markov model represents the clinical pathway using a set of mutually exclusive health states. At the beginning of every cycle, each person must be in exactly one state. Together, the states should cover all relevant situations that a person can experience within the model.

Our model contains seven health states:

Non-surgical treatment: the person receives either current or optimised non-surgical care.
Successful non-surgical treatment: symptoms are managed without knee replacement.
Total knee replacement surgery: the person undergoes primary knee replacement.
Successful total knee replacement: the person remains in the post-surgical state.
Revision surgery: the person undergoes surgery to revise the knee replacement.
Successful revision: the person remains in the post-revision state.
Death: the person has died and cannot move to another state.

The original model shown above also includes people who have not yet developed symptomatic knee osteoarthritis. We omit this state because our model begins with 400 people aged 50 years who already have knee osteoarthritis and are starting non-surgical treatment.

The model therefore begins with the following cohort distribution:

\[ \mathbf{m}_0 = (400,0,0,0,0,0,0). \]

The row vector $\mathbf{m}_0$ describes the number of people in each health state at the beginning of the model (cycle 0). Its seven positions correspond to the seven health states in the order defined above. All 400 people initially occupy NonSurgical; no one has yet entered another state.

The following hidden code creates the file containing the transition probabilities used later in the exercise. The parameters originate from the article (Table 3). files. You don’t have to run this code, you can import the input parameters later.

Code

transition.inputs <- data.frame(
  age = 50:89,
  mortality = c(
    0.001467, 0.001708, 0.002210, 0.002221, 0.002270,
    0.002567, 0.003225, 0.003103, 0.003706, 0.003826,
    0.004430, 0.004854, 0.005644, 0.006247, 0.007226,
    0.007951, 0.007912, 0.008952, 0.010508, 0.010098,
    0.012707, 0.013297, 0.015418, 0.016385, 0.019877,
    0.020639, 0.024836, 0.025394, 0.032061, 0.034113,
    0.038527, 0.044821, 0.051650, 0.056351, 0.066467,
    0.083008, 0.094663, 0.109233, 0.125404, 0.141392
  )
)

transition.inputs$p_tkr <- ifelse(
  transition.inputs$age <= 54,
  0.008,
  ifelse(
    transition.inputs$age <= 64,
    0.025,
    ifelse(
      transition.inputs$age <= 74,
      0.049,
      ifelse(
        transition.inputs$age <= 84,
        0.057,
        0.022
      )
    )
  )
)

transition.inputs$p_revision <- ifelse(
  transition.inputs$age <= 54,
  0.015625,
  ifelse(
    transition.inputs$age <= 64,
    0.011625,
    ifelse(
      transition.inputs$age <= 74,
      0.008125,
      ifelse(
        transition.inputs$age <= 84,
        0.005375,
        0.003500
      )
    )
  )
)

save(
  transition.inputs,
  file = "../data/markov-transition-inputs.rda"
)

rm(transition.inputs)

Exercise

1. Define the health states and starting cohort

Create a character vector named state.names containing the seven health states in the following order:

Code

state.names <- c(
  "NonSurgical",
  "SuccessfulNonSurgical",
  "TKR",
  "SuccessfulTKR",
  "Revision",
  "SuccessfulRevision",
  "Death"
)

Next, create the starting cohort:

Code

cohort.initial <- c(
  400, 0, 0, 0, 0, 0, 0
)

names(cohort.initial) <- state.names

Inspect the vector and use sum() to check that it contains 400 people:

Code

cohort.initial

sum(cohort.initial)

Consider why all participants begin in NonSurgical and why no one begins in TKR, Revision or Death.

2. Load and inspect the transition probabilities

The age-specific transition probabilities are stored in markov-transition-inputs.rda. Load the file using:

Code

load("../data/markov-transition-inputs.rda")

Remember that load() restores the object under the name with which it was saved. No assignment is required. After loading the file, the dataset is available as transition.inputs.

Inspect the dataset and its structure:

Code

transition.inputs

str(transition.inputs)

The dataset contains one row for every age from 50 to 89 and the following variables:

age: age at the beginning of the model cycle;
mortality: annual probability of death;
p_tkr: annual probability of moving from SuccessfulNonSurgical to TKR;
p_revision: annual probability of moving from SuccessfulTKR to Revision.

The model contains 40 transitions. The first occurs at age 50, and the final transition takes the cohort from age 89 to age 90.

Examine how the probabilities change with age. In particular, identify the ages at which p_tkr and p_revision change.

3. Create the transition-matrix function

A transition matrix contains the probabilities of moving between health states during one model cycle. Its rows represent the current state, while its columns represent the state in the following cycle.

In this model, $t$ denotes the model cycle and $a_t$ denotes the cohort’s age at the beginning of that cycle. Because the cohort begins at age 50 and one cycle represents one year:

\[ a_t = 50+t. \]

For example, cycle $t=0$ begins at age 50, cycle $t=1$ begins at age 51 and cycle $t=20$ begins at age 70.

If $\mathbf{m}_t$ is the cohort distribution at the beginning of cycle $t$, the distribution at the beginning of the next cycle is:

\[ \mathbf{m}_{t+1} = \mathbf{m}_t\mathbf{P}(a_t), \]

where:

$t$ is the model cycle;
$a_t$ is the cohort’s age at the beginning of cycle $t$;
$\mathbf{m}_t$ is the cohort vector (a row vector) at the beginning of cycle $t$;
$\mathbf{P}(a_t)$ is the transition matrix containing the probabilities for age $a_t$;
$\mathbf{m}_{t+1}$ is the cohort vector at the beginning of cycle $t+1$.

The first transition moves the cohort from age 50 to age 51:

\[ \mathbf{m}_1 = \mathbf{m}_0\mathbf{P}(50). \]

Using the transition probabilities for age 50, the transition matrix is (check with transition.inputs):

\[ \mathbf{P}(50) = \begin{pmatrix} 0 & 0.998533 & 0 & 0 & 0 & 0 & 0.001467 \\ 0 & 0.990533 & 0.008 & 0 & 0 & 0 & 0.001467 \\ 0 & 0 & 0 & 0.998533 & 0 & 0 & 0.001467 \\ 0 & 0 & 0 & 0.982908 & 0.015625 & 0 & 0.001467 \\ 0 & 0 & 0 & 0 & 0 & 0.998533 & 0.001467 \\ 0 & 0 & 0 & 0 & 0 & 0.998533 & 0.001467 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix}. \]

The rows and columns follow the same order:

NonSurgical;
SuccessfulNonSurgical;
TKR;
SuccessfulTKR;
Revision;
SuccessfulRevision;
Death.

Each cell is read from the state represented by its row to the state represented by its column. For example, the second row describes participants who are currently in SuccessfulNonSurgical:

\[ (0,\ 0.990533,\ 0.008,\ 0,\ 0,\ 0,\ 0.001467). \]

During the transition from age 50 to age 51, these participants have a probability of:

$0.990533$ of remaining in SuccessfulNonSurgical;
$0.008$ of moving to TKR;
$0.001467$ of moving to Death.

The remaining transitions have a probability of zero because they are not possible directly from SuccessfulNonSurgical.

The probabilities in each row must sum to 1. For the SuccessfulNonSurgical row:

\[ 0.990533+0.008+0.001467=1. \]

This ensures that every person in the state is assigned to exactly one state in the following cycle.

The transition probabilities are a function of age. This means, we now need a seperate transition matrix for each cycle. Of course, we don’t do this by hand. The following function creates the transition matrix for a specified age. Use the complete code as provided.

You are not expected to reproduce the function from memory. However, read through it and try to understand how the probabilities are selected and entered into the matrix. The figure above may help to understand it better.

Code

create_transition_matrix <- function(age) {
  
  age.row <- match(
    age,
    transition.inputs$age
  )
  
  p_death <- transition.inputs$mortality[
    age.row
  ]
  
  p_tkr <- transition.inputs$p_tkr[
    age.row
  ]
  
  p_revision <- transition.inputs$p_revision[
    age.row
  ]
  
  transition.matrix <- matrix(
    0,
    nrow = length(state.names),
    ncol = length(state.names),
    dimnames = list(
      from = state.names,
      to = state.names
    )
  )
  
  transition.matrix[
    "NonSurgical",
    "SuccessfulNonSurgical"
  ] <- 1 - p_death
  
  transition.matrix[
    "NonSurgical",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "SuccessfulNonSurgical",
    "SuccessfulNonSurgical"
  ] <- 1 - p_tkr - p_death
  
  transition.matrix[
    "SuccessfulNonSurgical",
    "TKR"
  ] <- p_tkr
  
  transition.matrix[
    "SuccessfulNonSurgical",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "TKR",
    "SuccessfulTKR"
  ] <- 1 - p_death
  
  transition.matrix[
    "TKR",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "SuccessfulTKR",
    "SuccessfulTKR"
  ] <- 1 - p_revision - p_death
  
  transition.matrix[
    "SuccessfulTKR",
    "Revision"
  ] <- p_revision
  
  transition.matrix[
    "SuccessfulTKR",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "Revision",
    "SuccessfulRevision"
  ] <- 1 - p_death
  
  transition.matrix[
    "Revision",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "SuccessfulRevision",
    "SuccessfulRevision"
  ] <- 1 - p_death
  
  transition.matrix[
    "SuccessfulRevision",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "Death",
    "Death"
  ] <- 1
  
  transition.matrix
}

Within the function, identify where the code:

finds the row for the specified age;
extracts the mortality, TKR and revision probabilities;
creates an empty $7\times7$ matrix;
enters the possible transitions between states;
calculates the probabilities of remaining in a state;
makes Death an absorbing state;
returns the completed matrix.

4. Inspect and validate the transition matrices

Use the function to create the transition matrix for age 50:

Code

transition.age50 <- create_transition_matrix(
  age = 50
)

Display the matrix rounded to four decimal places:

Code

round(transition.age50, 4)

Remember that:

each row represents a participant’s current state;
each column represents the participant’s state in the following cycle;
each cell contains the probability of the corresponding transition.

Use rowSums() to check that every row sums to 1:

Code

rowSums(transition.age50)

Finally, create another transition matrix for age 70:

Code

transition.age70 <- create_transition_matrix(
  age = 70
)

round(transition.age70, 4)

Compare the matrices for ages 50 and 70. Consider:

which transition probabilities differ;
why these probabilities differ;
how the changes affect the probability of remaining in the same state.

The two-year delay of total knee replacement under optimised care will be introduced in the next step, when we follow the cohort over repeated cycles.

1. Define the health states and starting cohort

We first create the vector containing the health-state names:

Code

state.names <- c(
  "NonSurgical",
  "SuccessfulNonSurgical",
  "TKR",
  "SuccessfulTKR",
  "Revision",
  "SuccessfulRevision",
  "Death"
)

The order is important because it will also be used for the cohort vectors and the rows and columns of every transition matrix.

We then create the initial cohort:

Code

cohort.initial <- c(
  400, 0, 0, 0, 0, 0, 0
)

names(cohort.initial) <- state.names

cohort.initial

          NonSurgical SuccessfulNonSurgical                   TKR 
                  400                     0                     0 
        SuccessfulTKR              Revision    SuccessfulRevision 
                    0                     0                     0 
                Death 
                    0

All 400 people begin in NonSurgical because the cohort consists of people with knee osteoarthritis who are starting non-surgical treatment. No one has yet undergone knee replacement or revision surgery.

We verify the size of the cohort:

Code

sum(cohort.initial)

[1] 400

The sum is 400, as expected.

2. Load and inspect the transition probabilities

We load the prepared model inputs:

Code

load("../data/markov-transition-inputs.rda")

The file restores the object named transition.inputs. It already contains the age-specific mortality, total knee replacement and revision probabilities, so we do not need to calculate these inputs ourselves.

We inspect the structure of the dataset:

Code

str(transition.inputs)

'data.frame':   40 obs. of  4 variables:
 $ age       : int  50 51 52 53 54 55 56 57 58 59 ...
 $ mortality : num  0.00147 0.00171 0.00221 0.00222 0.00227 ...
 $ p_tkr     : num  0.008 0.008 0.008 0.008 0.008 0.025 0.025 0.025 0.025 0.025 ...
 $ p_revision: num  0.0156 0.0156 0.0156 0.0156 0.0156 ...

The dataset contains 40 observations and four variables:

age;
mortality;
p_tkr;
p_revision.

There should be no missing values because every model cycle requires a complete set of transition probabilities.

Mortality changes with every year of age. The probabilities of total knee replacement and revision surgery remain constant within defined age groups.

We inspect the rows around the boundaries of these age groups:

Code

transition.inputs[
  c(1, 5, 6, 15, 16, 25, 26, 35, 36, 40),
]

   age mortality p_tkr p_revision
1   50  0.001467 0.008   0.015625
5   54  0.002270 0.008   0.015625
6   55  0.002567 0.025   0.011625
15  64  0.007226 0.025   0.011625
16  65  0.007951 0.049   0.008125
25  74  0.019877 0.049   0.008125
26  75  0.020639 0.057   0.005375
35  84  0.066467 0.057   0.005375
36  85  0.083008 0.022   0.003500
40  89  0.141392 0.022   0.003500

The probability of total knee replacement increases across the age groups until age 84 and then decreases in the oldest age group. The annual probability of revision surgery decreases across the older age groups.

3. Create the transition-matrix function

We use the function provided in the task:

Code

create_transition_matrix <- function(age) {
  
  age.row <- match(
    age,
    transition.inputs$age
  )
  
  p_death <- transition.inputs$mortality[
    age.row
  ]
  
  p_tkr <- transition.inputs$p_tkr[
    age.row
  ]
  
  p_revision <- transition.inputs$p_revision[
    age.row
  ]
  
  transition.matrix <- matrix(
    0,
    nrow = length(state.names),
    ncol = length(state.names),
    dimnames = list(
      from = state.names,
      to = state.names
    )
  )
  
  transition.matrix[
    "NonSurgical",
    "SuccessfulNonSurgical"
  ] <- 1 - p_death
  
  transition.matrix[
    "NonSurgical",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "SuccessfulNonSurgical",
    "SuccessfulNonSurgical"
  ] <- 1 - p_tkr - p_death
  
  transition.matrix[
    "SuccessfulNonSurgical",
    "TKR"
  ] <- p_tkr
  
  transition.matrix[
    "SuccessfulNonSurgical",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "TKR",
    "SuccessfulTKR"
  ] <- 1 - p_death
  
  transition.matrix[
    "TKR",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "SuccessfulTKR",
    "SuccessfulTKR"
  ] <- 1 - p_revision - p_death
  
  transition.matrix[
    "SuccessfulTKR",
    "Revision"
  ] <- p_revision
  
  transition.matrix[
    "SuccessfulTKR",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "Revision",
    "SuccessfulRevision"
  ] <- 1 - p_death
  
  transition.matrix[
    "Revision",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "SuccessfulRevision",
    "SuccessfulRevision"
  ] <- 1 - p_death
  
  transition.matrix[
    "SuccessfulRevision",
    "Death"
  ] <- p_death
  
  transition.matrix[
    "Death",
    "Death"
  ] <- 1
  
  transition.matrix
}

The object age.row contains the position of the specified age in the imported dataset. For example:

Code

match(
  70,
  transition.inputs$age
)

[1] 21

This returns the row containing the transition probabilities for age 70. The function then extracts all three probabilities from this row.

The code creates a $7\times7$ matrix containing zeros. The state names are assigned to its rows and columns before the transition probabilities are entered.

For example:

Code

transition.matrix[
  "SuccessfulNonSurgical",
  "TKR"
] <- p_tkr

This places the age-specific probability of total knee replacement in the row representing the current state and the column representing the next state.

The probability of remaining in SuccessfulNonSurgical is therefore calculated as:

\[ P(\text{remain in SuccessfulNonSurgical}) = 1-p_{\text{TKR}}-p_{\text{death}}. \]

The temporary states NonSurgical, TKR and Revision do not have self-transitions. Participants who survive the cycle move to the corresponding successful state.

Finally, the following assignment makes Death an absorbing state:

Code

transition.matrix[
  "Death",
  "Death"
] <- 1

This means that a person who has entered Death remains there in every subsequent cycle:

\[ P(\text{Death}\rightarrow\text{Death})=1. \]

All probabilities of moving from Death to another state remain zero.

4. Inspect and validate the transition matrices

We first create the transition matrix for age 50:

Code

transition.age50 <- create_transition_matrix(
  age = 50
)

round(transition.age50, 4)

                       to
from                    NonSurgical SuccessfulNonSurgical   TKR SuccessfulTKR
  NonSurgical                     0                0.9985 0.000        0.0000
  SuccessfulNonSurgical           0                0.9905 0.008        0.0000
  TKR                             0                0.0000 0.000        0.9985
  SuccessfulTKR                   0                0.0000 0.000        0.9829
  Revision                        0                0.0000 0.000        0.0000
  SuccessfulRevision              0                0.0000 0.000        0.0000
  Death                           0                0.0000 0.000        0.0000
                       to
from                    Revision SuccessfulRevision  Death
  NonSurgical             0.0000             0.0000 0.0015
  SuccessfulNonSurgical   0.0000             0.0000 0.0015
  TKR                     0.0000             0.0000 0.0015
  SuccessfulTKR           0.0156             0.0000 0.0015
  Revision                0.0000             0.9985 0.0015
  SuccessfulRevision      0.0000             0.9985 0.0015
  Death                   0.0000             0.0000 1.0000

Each row describes what can happen to participants who begin the cycle in that state.

For example, participants in SuccessfulNonSurgical can:

remain in SuccessfulNonSurgical;
move to TKR;
move to Death.

All other transitions from this state have a probability of zero.

We check the row sums:

Code

rowSums(transition.age50)

          NonSurgical SuccessfulNonSurgical                   TKR 
                    1                     1                     1 
        SuccessfulTKR              Revision    SuccessfulRevision 
                    1                     1                     1 
                Death 
                    1

Every row sums to 1. This confirms that all possible outcomes during the cycle have been included.

We then create the transition matrix for age 70:

Code

transition.age70 <- create_transition_matrix(
  age = 70
)

round(transition.age70, 4)

                       to
from                    NonSurgical SuccessfulNonSurgical   TKR SuccessfulTKR
  NonSurgical                     0                0.9873 0.000        0.0000
  SuccessfulNonSurgical           0                0.9383 0.049        0.0000
  TKR                             0                0.0000 0.000        0.9873
  SuccessfulTKR                   0                0.0000 0.000        0.9792
  Revision                        0                0.0000 0.000        0.0000
  SuccessfulRevision              0                0.0000 0.000        0.0000
  Death                           0                0.0000 0.000        0.0000
                       to
from                    Revision SuccessfulRevision  Death
  NonSurgical             0.0000             0.0000 0.0127
  SuccessfulNonSurgical   0.0000             0.0000 0.0127
  TKR                     0.0000             0.0000 0.0127
  SuccessfulTKR           0.0081             0.0000 0.0127
  Revision                0.0000             0.9873 0.0127
  SuccessfulRevision      0.0000             0.9873 0.0127
  Death                   0.0000             0.0000 1.0000

We compare the transitions from successful non-surgical treatment:

Code

transition.age50[
  "SuccessfulNonSurgical",
]

          NonSurgical SuccessfulNonSurgical                   TKR 
             0.000000              0.990533              0.008000 
        SuccessfulTKR              Revision    SuccessfulRevision 
             0.000000              0.000000              0.000000 
                Death 
             0.001467

Code

transition.age70[
  "SuccessfulNonSurgical",
]

          NonSurgical SuccessfulNonSurgical                   TKR 
             0.000000              0.938293              0.049000 
        SuccessfulTKR              Revision    SuccessfulRevision 
             0.000000              0.000000              0.000000 
                Death 
             0.012707

At age 50, the annual probability of total knee replacement is 0.008. At age 70, it is 0.049. Mortality is also higher at age 70.

Consequently, the probability of remaining in SuccessfulNonSurgical is lower at age 70:

\[ P(\text{remain in SuccessfulNonSurgical}) = 1-p_{\text{TKR}}-p_{\text{death}}. \]

We also compare the transitions following successful knee replacement:

Code

transition.age50[
  "SuccessfulTKR",
]

          NonSurgical SuccessfulNonSurgical                   TKR 
             0.000000              0.000000              0.000000 
        SuccessfulTKR              Revision    SuccessfulRevision 
             0.982908              0.015625              0.000000 
                Death 
             0.001467

Code

transition.age70[
  "SuccessfulTKR",
]

          NonSurgical SuccessfulNonSurgical                   TKR 
             0.000000              0.000000              0.000000 
        SuccessfulTKR              Revision    SuccessfulRevision 
             0.979168              0.008125              0.000000 
                Death 
             0.012707

The age-specific revision and mortality probabilities differ. A new transition matrix must therefore be created during every model cycle as the cohort grows older. So there is indeed a reason for the nurd-stuff :-).

Competing transitions

Total knee replacement, revision surgery and death are treated as competing transitions within a cycle. The probability of remaining in the current state is calculated after subtracting all probabilities of leaving that state.

For example:

\[ P(\text{remain in SuccessfulTKR}) = 1-p_{\text{revision}}-p_{\text{death}}. \]

This ensures that every row sums to 1 and prevents the same person from making two different transitions during the same cycle.

3.2 2. Follow the cohort over repeated cycles

A transition matrix describes what may happen during one cycle. To project the cohort over the complete 40-year time horizon, we apply the appropriate age-specific transition matrix repeatedly.

The initial cohort distribution is:

\[ \mathbf{m}_0 = (400,0,0,0,0,0,0). \]

The first transition uses the matrix for age 50:

\[ \mathbf{m}_1 = \mathbf{m}_0\mathbf{P}(50). \]

The resulting row vector $\mathbf{m}_1$ describes the cohort at age 51. It is then multiplied by the matrix for age 51:

\[ \mathbf{m}_2 = \mathbf{m}_1\mathbf{P}(51). \]

More generally:

\[ \mathbf{m}_{t+1} = \mathbf{m}_t\mathbf{P}(a_t), \qquad a_t=50+t, \]

where:

$t$ is the model cycle;
$a_t$ is the cohort’s age at the beginning of cycle $t$;
$\mathbf{m}_t$ is the cohort distribution at the beginning of cycle $t$;
$\mathbf{P}(a_t)$ is the transition matrix for age $a_t$;
$\mathbf{m}_{t+1}$ is the cohort distribution at the beginning of the next cycle.

Repeating this calculation produces a cohort trace. A cohort trace is a matrix in which:

each row represents one model cycle;
each column represents one health state;
each cell contains the expected number of people occupying that state.

Because transition probabilities are applied to a hypothetical cohort, the state counts do not need to be whole numbers. For example, the model may predict that 3.4 people undergo total knee replacement during a cycle. This does not describe 0.4 of an actual person; it represents the expected number in a cohort of 400.

The trace contains 41 rows:

row 1 represents cycle 0 at age 50, before any transition;
row 2 represents cycle 1 at age 51;
the final row represents cycle 40 at age 90.

There are only 40 transitions because the starting distribution is recorded before the first transition.

Representing the two-year delay

Under current care, total knee replacement becomes possible once participants have completed the initial cycle of non-surgical treatment and entered SuccessfulNonSurgical.

Under optimised care, we assume that total knee replacement is delayed by two additional years. We therefore prevent transitions from SuccessfulNonSurgical to TKR during the first two cycles in which knee replacement would otherwise be possible:

during the transition from age 51 to age 52;
during the transition from age 52 to age 53.

During these two transitions, the probability of moving to TKR is set to zero. The corresponding probability is added to the probability of remaining in SuccessfulNonSurgical, ensuring that the row still sums to 1. From the transition beginning at age 53 onwards, total knee replacement is permitted again.

Simplifying the delay

The original model represented delayed knee replacement using additional tunnel states. Tunnel states retain information about how long people have spent in a particular part of the model.

For this teaching exercise, everyone enters non-surgical treatment at the same age. We can therefore represent the two-year delay by temporarily preventing the transition to TKR. This is easier to implement, but it is not identical to shifting every individual knee replacement exactly two years into the future.

The simplification allows us to focus on the main cohort calculations without inviting tunnel states into the exercise.

For both strategies, the complete cohort remains in the model. As people die, they accumulate in the absorbing Death state. Consequently, the sum across all seven states should remain equal to 400 in every cycle:

\[ \sum_{j=1}^{7}m_{tj}=400. \]

This provides a useful check that no members of the cohort have quietly disappeared during the calculations.

Exercise

1. Prepare the cohort traces

The model covers 40 annual cycles, beginning at age 50 and ending at age 90. Define the model settings:

Code

start.age <- 50
n.cycles <- 40

model.ages <- start.age:(start.age + n.cycles)

Next, create an empty matrix for each strategy:

Code

cohort.current <- matrix(
  0,
  nrow = n.cycles + 1,
  ncol = length(state.names),
  dimnames = list(
    cycle = 0:n.cycles,
    state = state.names
  )
)

cohort.optimised <- matrix(
  0,
  nrow = n.cycles + 1,
  ncol = length(state.names),
  dimnames = list(
    cycle = 0:n.cycles,
    state = state.names
  )
)

Why do these matrices contain 41 rows even though the model has only 40 cycles?

Enter the initial cohort distribution into the first row of both matrices:

Code

cohort.current[1, ] <- cohort.initial
cohort.optimised[1, ] <- cohort.initial

Both strategies begin with the same cohort. The strategies only differ in what happens after the model starts. You can look at the objects cohort.current and cohort.optimised the see matrices.

2. Follow the cohort under current care

The following loop applies the appropriate age-specific transition matrix during each cycle:

Code

for (cycle in 1:n.cycles) {
  
  current.age <- start.age + cycle - 1
  
  transition.current <-
    create_transition_matrix(
      age = current.age
    )
  
  cohort.current[cycle + 1, ] <-
    cohort.current[cycle, ] %*%
    transition.current
}

In the code above, %*% is the operator for matrix multiplication. So the last bit

Code

cohort.current[cycle + 1, ] <-
    cohort.current[cycle, ] %*%
    transition.current

is equivalent with

\[ \mathbf{m}_{t+1} = \mathbf{m}_t\mathbf{P}(a_t). \]

Run the code and inspect the first six rows:

Code

round(
  cohort.current[1:6, ],
  2
)

In each cycle, it multiplies the current cohort vector by the transition matrix for the relevant age.

Use rowSums() to check that the total cohort remains equal to 400 throughout the model:

Code

round(
  rowSums(cohort.current),
  6
)

3. Follow the cohort under optimised care

The optimised strategy uses the same transition probabilities, except that total knee replacement is prevented during the transitions beginning at ages 51 and 52.

During these two cycles, the probability of moving to TKR must be:

obtained from the transition matrix;
added to the probability of remaining in SuccessfulNonSurgical;
replaced by zero in the TKR column.

The required modification is:

Code

# 1: 
delayed.probability <- transition.optimised[ 
  "SuccessfulNonSurgical",
  "TKR"
]

# 2:
transition.optimised[ 
  "SuccessfulNonSurgical",
  "SuccessfulNonSurgical"
] <- transition.optimised[
  "SuccessfulNonSurgical",
  "SuccessfulNonSurgical"
] + delayed.probability

# 3:
transition.optimised[
  "SuccessfulNonSurgical",
  "TKR"
] <- 0

Place this modification inside an if statement so that it is applied only at ages 51 and 52:

Code

for (cycle in 1:n.cycles) {
  
  current.age <- start.age + cycle - 1
  
  transition.optimised <-
    create_transition_matrix(
      age = current.age
    )
  
  if (current.age %in% c(51, 52)) {
    
    # Add modification to temporarily prevent the transition to TKR
    
  }
  
  cohort.optimised[cycle + 1, ] <-
    cohort.optimised[cycle, ] %*%
    transition.optimised
}

Run the loop from above. Then, use rowSums() to confirm that the optimised-care trace also contains 400 people in every row.

4. Add cycle and age information

The row names show the model cycles, but adding explicit cycle and age variables makes the cohort traces easier to inspect.

Convert the two matrices into data frames:

Code

trace.current <- data.frame(
  cycle = 0:n.cycles,
  age = model.ages,
  cohort.current,
  check.names = FALSE
)

trace.optimised <- data.frame(
  cycle = 0:n.cycles,
  age = model.ages,
  cohort.optimised,
  check.names = FALSE
)

Inspect the both traces using head() and tail() as they are data frames now.

5. Compare the cohort traces

First, inspect the beginning of the two cohort traces directly. This makes it easier to identify when total knee replacement first occurs under each strategy. Use:

Code

trace.current[
  1:8,
  c(
    "cycle",
    "age",
    "SuccessfulNonSurgical",
    "TKR"
  )
]

trace.optimised[
  1:8,
  c(
    "cycle",
    "age",
    "SuccessfulNonSurgical",
    "TKR"
  )
]

Compare the values and identify:

the first age at which participants enter TKR under current care;
the first age at which participants enter TKR under optimised care;
what happens to the number of people in SuccessfulNonSurgical during the delay.

Next, compare current and optimised care for the four persistent health states. Create one plot for each state. Each plot should contain:

a black line for current care;
a blue line for optimised care.

Use the following code:

Code

states.to.plot <- c(
  "SuccessfulNonSurgical",
  "SuccessfulTKR",
  "SuccessfulRevision",
  "Death"
)

old.par <- par(
  no.readonly = TRUE
)

par(
  mfrow = c(2, 2),
  mar = c(4, 4, 3, 1)
)

for (state in states.to.plot) {
  
  y.range <- range(
    c(
      trace.current[[state]],
      trace.optimised[[state]]
    )
  )
  
  plot(
    trace.current$age,
    trace.current[[state]],
    type = "l",
    lwd = 2,
    col = "black",
    ylim = y.range,
    main = state,
    xlab = "Age",
    ylab = "Number of people"
  )
  
  lines(
    trace.optimised$age,
    trace.optimised[[state]],
    lwd = 2,
    col = "blue"
  )
  
  legend(
    "topright",
    legend = c(
      "Current",
      "Optimised"
    ),
    col = c(
      "black",
      "blue"
    ),
    lty = 1,
    lwd = 2,
    cex = 0.7,
    bty = "n"
  )
}

par(old.par)

The expression trace.current[[state]] selects the column whose name is currently stored in state. The loop repeats the plotting commands for each of the four states listed in states.to.plot.

Use the plots to examine:

when differences between the strategies first appear;
how the delay affects SuccessfulNonSurgical and SuccessfulTKR;
whether the delay reduces the number of people reaching SuccessfulRevision;
whether the number of deaths differs between the strategies.

Why are the surgical states not plotted?

TKR and Revision are temporary states. They show the expected number of procedures occurring during a particular cycle rather than the number of people remaining in a long-term health state.

These relatively small annual numbers are easier to compare directly in the cohort traces. NonSurgical is also omitted because everyone leaves this initial state after the first cycle.

1. Prepare the cohort traces

We first define the starting age, number of cycles and ages represented in the cohort traces:

Code

start.age <- 50
n.cycles <- 40

model.ages <- start.age:(start.age + n.cycles)

The vector model.ages contains 41 ages, from 50 to 90. It includes the starting distribution at age 50 and the distribution after the final transition at age 90.

We create an empty matrix for each strategy:

Code

cohort.current <- matrix(
  0,
  nrow = n.cycles + 1,
  ncol = length(state.names),
  dimnames = list(
    cycle = 0:n.cycles,
    state = state.names
  )
)

cohort.optimised <- matrix(
  0,
  nrow = n.cycles + 1,
  ncol = length(state.names),
  dimnames = list(
    cycle = 0:n.cycles,
    state = state.names
  )
)

Each matrix has:

41 rows representing cycle 0 to cycle 40;
seven columns representing the health states.

The model has 40 transitions but 41 recorded cohort distributions because the initial distribution is recorded before the first transition.

We enter the starting cohort into the first row of both matrices:

Code

cohort.current[1, ] <- cohort.initial
cohort.optimised[1, ] <- cohort.initial

cohort.current[1, ]

          NonSurgical SuccessfulNonSurgical                   TKR 
                  400                     0                     0 
        SuccessfulTKR              Revision    SuccessfulRevision 
                    0                     0                     0 
                Death 
                    0

Code

cohort.optimised[1, ]

          NonSurgical SuccessfulNonSurgical                   TKR 
                  400                     0                     0 
        SuccessfulTKR              Revision    SuccessfulRevision 
                    0                     0                     0 
                Death 
                    0

Both strategies begin with 400 people in NonSurgical.

2. Follow the cohort under current care

We now apply the appropriate transition matrix repeatedly:

Code

for (cycle in 1:n.cycles) {
  
  current.age <- start.age + cycle - 1
  
  transition.current <-
    create_transition_matrix(
      age = current.age
    )
  
  cohort.current[cycle + 1, ] <-
    cohort.current[cycle, ] %*%
    transition.current
}

In the first repetition, cycle equals 1 and current.age is:

\[ 50+1-1=50. \]

The calculation performed is therefore:

\[ \mathbf{m}_1 = \mathbf{m}_0\mathbf{P}(50). \]

The result is stored in row 2 because row 1 contains the starting distribution.

In the second repetition, the model uses the transition matrix for age 51:

\[ \mathbf{m}_2 = \mathbf{m}_1\mathbf{P}(51). \]

This continues until the transition from age 89 to age 90 has been completed.

We inspect the first six rows:

Code

round(
  cohort.current[1:6, ],
  2
)

     state
cycle NonSurgical SuccessfulNonSurgical  TKR SuccessfulTKR Revision
    0         400                  0.00 0.00          0.00     0.00
    1           0                399.41 0.00          0.00     0.00
    2           0                395.54 3.20          0.00     0.00
    3           0                391.50 3.16          3.19     0.00
    4           0                387.50 3.13          6.29     0.05
    5           0                383.52 3.10          9.30     0.10
     state
cycle SuccessfulRevision Death
    0               0.00  0.00
    1               0.00  0.59
    2               0.00  1.27
    3               0.00  2.15
    4               0.00  3.03
    5               0.05  3.93

After the first transition, almost all surviving participants have moved from NonSurgical to SuccessfulNonSurgical. During subsequent cycles, some remain in non-surgical care, some undergo total knee replacement and some die.

The state counts are not necessarily whole numbers. They represent expected numbers in the hypothetical cohort rather than counts of individually simulated people.

We check that the total cohort remains equal to 400:

Code

round(
  rowSums(cohort.current),
  6
)

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19 
400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 
 20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39 
400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 400 
 40 
400

Every row sums to 400. People who die remain in the Death state, so they are still included in the cohort trace.

A more concise check is:

Code

all(
  abs(rowSums(cohort.current) - 400) < 0.000001
)

[1] TRUE

The result is TRUE, confirming that no one has been lost from the model. This is always a very important check regarding technical validity.

3. Follow the cohort under optimised care

We construct the optimised-care trace using the same process. During the transitions beginning at ages 51 and 52, however, we temporarily prevent total knee replacement:

Code

for (cycle in 1:n.cycles) {
  
  current.age <- start.age + cycle - 1
  
  transition.optimised <-
    create_transition_matrix(
      age = current.age
    )
  
  if (current.age %in% c(51, 52)) {
    
    delayed.probability <- transition.optimised[
      "SuccessfulNonSurgical",
      "TKR"
    ]
    
    transition.optimised[
      "SuccessfulNonSurgical",
      "SuccessfulNonSurgical"
    ] <- transition.optimised[
      "SuccessfulNonSurgical",
      "SuccessfulNonSurgical"
    ] + delayed.probability
    
    transition.optimised[
      "SuccessfulNonSurgical",
      "TKR"
    ] <- 0
  }
  
  cohort.optimised[cycle + 1, ] <-
    cohort.optimised[cycle, ] %*%
    transition.optimised
}

The transition probability that would normally move people to TKR is added to the probability of remaining in SuccessfulNonSurgical. This keeps the row sum equal to 1.

For example, at age 51:

\[ P(\text{remain in SuccessfulNonSurgical}) = 1-p_{\text{death}} \]

because $p_{\text{TKR}}$ is temporarily set to zero.

We verify that the cohort is preserved:

Code

all(
  abs(rowSums(cohort.optimised) - 400) < 0.000001
)

[1] TRUE

The result is again TRUE 🍺 .

Why ages 51 and 52?

During the first transition, from age 50 to age 51, everyone begins in NonSurgical. The model does not permit a direct transition from NonSurgical to TKR.

Under current care, the first opportunity for TKR therefore occurs during the transition from age 51 to age 52. Preventing TKR during the transitions beginning at ages 51 and 52 delays the first possible knee replacements by two years. Under optimised care, TKR first becomes possible during the transition from age 53 to age 54.

4. Add cycle and age information

We convert the cohort matrices into data frames and add the corresponding cycle and age:

Code

trace.current <- data.frame(
  cycle = 0:n.cycles,
  age = model.ages,
  cohort.current,
  check.names = FALSE
)

trace.optimised <- data.frame(
  cycle = 0:n.cycles,
  age = model.ages,
  cohort.optimised,
  check.names = FALSE
)

The argument check.names = FALSE preserves the health-state names exactly as defined.

We inspect the beginning of the traces:

Code

head(trace.current)

  cycle age NonSurgical SuccessfulNonSurgical      TKR SuccessfulTKR   Revision
0     0  50         400                0.0000 0.000000      0.000000 0.00000000
1     1  51           0              399.4132 0.000000      0.000000 0.00000000
2     2  52           0              395.5357 3.195306      0.000000 0.00000000
3     3  53           0              391.4973 3.164286      3.188244 0.00000000
4     4  54           0              387.4958 3.131978      6.288604 0.04981631
5     5  55           0              383.5162 3.099966      9.300938 0.09825944
  SuccessfulRevision    Death
0         0.00000000 0.000000
1         0.00000000 0.586800
2         0.00000000 1.268998
3         0.00000000 2.150193
4         0.00000000 3.033818
5         0.04970323 3.934931

Code

head(trace.optimised)

  cycle age NonSurgical SuccessfulNonSurgical      TKR SuccessfulTKR Revision
0     0  50         400                0.0000 0.000000      0.000000        0
1     1  51           0              399.4132 0.000000      0.000000        0
2     2  52           0              398.7310 0.000000      0.000000        0
3     3  53           0              397.8498 0.000000      0.000000        0
4     4  54           0              393.7834 3.182798      0.000000        0
5     5  55           0              389.7392 3.150267      3.175574        0
  SuccessfulRevision    Death
0                  0 0.000000
1                  0 0.586800
2                  0 1.268998
3                  0 2.150193
4                  0 3.033818
5                  0 3.934931

We can also inspect the final cycles:

Code

tail(trace.current)

   cycle age NonSurgical SuccessfulNonSurgical       TKR SuccessfulTKR
35    35  85           0              54.67632 3.5555421     148.68119
36    36  86           0              48.93486 1.2028789     139.07948
37    37  87           0              43.22598 1.0765670     126.51603
38    38  88           0              37.55330 0.9509715     113.21247
39    39  89           0              32.01780 0.8261726      99.45065
40    40  90           0              26.78634 0.7043915      85.75040
    Revision SuccessfulRevision    Death
35 0.8393397           15.59296 176.6547
36 0.5203842           15.06829 195.1941
37 0.4867782           14.11300 214.5816
38 0.4428061           13.00500 234.8354
39 0.3962437           11.76140 255.5477
40 0.3480773           10.43865 275.9721

Code

tail(trace.optimised)

   cycle age NonSurgical SuccessfulNonSurgical       TKR SuccessfulTKR
35    35  85           0              55.56351 3.6132352     148.35437
36    36  86           0              49.72889 1.2223971     138.83384
37    37  87           0              43.92737 1.0940357     126.31217
38    38  88           0              38.16265 0.9664022     113.04715
39    39  89           0              32.53732 0.8395783      99.32013
40    40  90           0              27.22099 0.7158211      85.65031
    Revision SuccessfulRevision    Death
35 0.8370952           14.97714 176.6547
36 0.5192403           14.50153 195.1941
37 0.4859184           13.59886 214.5816
38 0.4420926           12.54625 234.8354
39 0.3956650           11.35956 255.5477
40 0.3476205           10.09313 275.9721

The Death column increases over time because mortality becomes more likely as the cohort ages. The number of living people therefore decreases, although the total across all seven states remains 400.

Note: the end of one cycle is simultaneously the beginning of the next cycle. Therefore, age 51 is both:

the end age of cycle 1
the starting age of cycle 2.

5. Compare the cohort traces

We first inspect the early model cycles:

Code

trace.current[
  1:8,
  c(
    "cycle",
    "age",
    "SuccessfulNonSurgical",
    "TKR"
  )
]

  cycle age SuccessfulNonSurgical      TKR
0     0  50                0.0000 0.000000
1     1  51              399.4132 0.000000
2     2  52              395.5357 3.195306
3     3  53              391.4973 3.164286
4     4  54              387.4958 3.131978
5     5  55              383.5162 3.099966
6     6  56              372.9438 9.587905
7     7  57              362.4175 9.323595

Code

trace.optimised[
  1:8,
  c(
    "cycle",
    "age",
    "SuccessfulNonSurgical",
    "TKR"
  )
]

  cycle age SuccessfulNonSurgical      TKR
0     0  50                0.0000 0.000000
1     1  51              399.4132 0.000000
2     2  52              398.7310 0.000000
3     3  53              397.8498 0.000000
4     4  54              393.7834 3.182798
5     5  55              389.7392 3.150267
6     6  56              378.9953 9.743481
7     7  57              368.2981 9.474882

Under current care, participants first enter TKR at age 52. This reflects transitions occurring between ages 51 and 52.

Under optimised care:

no participants enter TKR at age 52;
no participants enter TKR at age 53;
participants first enter TKR at age 54.

During the two-year delay, more people remain in SuccessfulNonSurgical under optimised care.

We now compare the four selected states graphically:

Code

states.to.plot <- c(
  "SuccessfulNonSurgical",
  "SuccessfulTKR",
  "SuccessfulRevision",
  "Death"
)

old.par <- par(
  no.readonly = TRUE
)

par(
  mfrow = c(2, 2),
  mar = c(4, 4, 3, 1)
)

for (state in states.to.plot) {
  
  y.range <- range(
    c(
      trace.current[[state]],
      trace.optimised[[state]]
    )
  )
  
  plot(
    trace.current$age,
    trace.current[[state]],
    type = "l",
    lwd = 2,
    col = "black",
    ylim = y.range,
    main = state,
    xlab = "Age",
    ylab = "Number of people"
  )
  
  lines(
    trace.optimised$age,
    trace.optimised[[state]],
    lwd = 2,
    col = "blue"
  )
  
  legend(
    "topright",
    legend = c(
      "Current",
      "Optimised"
    ),
    col = c(
      "black",
      "blue"
    ),
    lty = 1,
    lwd = 2,
    cex = 0.7,
    bty = "n"
  )
}

Code

par(old.par)

The differences first appear in SuccessfulNonSurgical. During the delay, more people remain in this state under optimised care because they cannot yet move to TKR.

The delay also changes the number of people in SuccessfulTKR. Because surgery occurs later, participants spend less time in the post-TKR state over the modelled period.

Fewer or later knee replacements also reduce the opportunity for subsequent revision surgery. The optimised strategy therefore results in fewer people reaching SuccessfulRevision.

The lines for Death overlap. This is expected because the model applies the same age-specific mortality probability to every living state under both strategies. Delaying TKR changes where living participants are located, but it does not change their mortality risk in this simplified model.

Why are the surgical states not plotted?

TKR and Revision are temporary states. They show the expected number of procedures occurring during a particular cycle rather than the number of people remaining in a long-term health state.

These relatively small annual numbers are easier to compare directly in the cohort traces. NonSurgical is also omitted because everyone leaves this initial state after the first cycle.

3.3 3. Assign costs and utilities to the health states

The cohort traces show how many people occupy each health state during every model cycle. To translate these clinical pathways into economic outcomes, we assign a cost and a utility value to each state.

A state cost represents the healthcare resources associated with occupying a health state for one cycle. For example:

NonSurgical includes the cost of the initial non-surgical treatment;
TKR includes the cost of total knee replacement surgery;
SuccessfulTKR includes the annual healthcare costs incurred after surgery;
Revision includes the cost of revision surgery;
Death has no further healthcare costs in this simplified model.

A state utility represents health-related quality of life while occupying that state. Utility values are measured on a scale on which:

$1$ represents full health;
$0$ represents a health state equivalent to death;
values between $0$ and $1$ represent different levels of health-related quality of life.

Because each model cycle lasts one year, occupying a state with utility $u$ for one complete cycle generates:

\[ QALY=u\times1=u. \]

For example, one year in a state with a utility of $0.75$ generates:

\[ 0.75\times1=0.75\text{ QALYs}. \]

The costs and utilities used in this exercise are based on the published evaluation by Vetsch et al. (2023) (Tables 1 & 2). Costs are expressed in 2019 Swiss francs.

Health state	Current-care cost	Optimised-care cost	Current-care utility	Optimised-care utility
`NonSurgical`	CHF 222	CHF 1,209	0.658	0.658
`SuccessfulNonSurgical`	CHF 200	CHF 200	0.749	0.783
`TKR`	CHF 18,326	CHF 18,326	0.658	0.658
`SuccessfulTKR`	CHF 500	CHF 500	0.878	0.878
`Revision`	CHF 28,776	CHF 28,776	0.658	0.658
`SuccessfulRevision`	CHF 500	CHF 500	0.760	0.760
`Death`	CHF 0	CHF 0	0	0

The strategies differ in two state values:

optimised non-surgical treatment has a higher initial cost than current non-surgical treatment;
SuccessfulNonSurgical has a higher utility under optimised care because the treatment is assumed to produce a larger improvement in health-related quality of life.

All other state costs and utilities are identical between the strategies. Nevertheless, the total costs and QALYs may differ because the strategies distribute the cohort differently across the states over time.

For example, delaying TKR may:

increase the time spent in SuccessfulNonSurgical;
reduce the time spent in SuccessfulTKR;
reduce the number of revision procedures;
avoid some expensive surgical costs.

The state values can be stored as cost and utility vectors. Their order must match the order of the columns in the cohort traces:

\[ \mathbf{c} = (c_1,c_2,\ldots,c_7) \]

and:

\[ \mathbf{u} = (u_1,u_2,\ldots,u_7), \]

where:

$c_j$ is the annual cost assigned to health state $j$;
$u_j$ is the utility assigned to health state $j$.

For a cohort distribution $\mathbf{m}_t$, the costs generated during cycle $t$ are:

\[ C_t = \mathbf{m}_t\mathbf{c}, \]

and the QALYs generated during that cycle are:

\[ Q_t = \mathbf{m}_t\mathbf{u}. \]

These calculations multiply the number of people in each state by the corresponding state cost or utility and then sum across all states.

Timing assumption

We assume that the cohort distribution at the beginning of a cycle determines the costs and utilities accrued during that cycle. Each person is assigned the full annual cost and utility of their current state.

The model does not apply a half-cycle correction. This is a teaching simplification: real models may account more precisely for transitions occurring partway through a cycle.

At this stage, we will define the state costs and utilities and calculate the undiscounted values generated during each cycle. Discounting and accumulation over the complete time horizon follow in the next step.

Exercise

1. Create the state-cost vectors

Create a (named) vector containing the annual state costs for current care. The order of the values must match the order of state.names. Use:

Code

costs.current <- c(
  NonSurgical = 222,
  SuccessfulNonSurgical = 200,
  TKR = 18326,
  SuccessfulTKR = 500,
  Revision = 28776,
  SuccessfulRevision = 500,
  Death = 0
)

Next, create a vector named costs.optimised for optimised care. Use the same costs, except that the cost of NonSurgical is CHF 1,209.

You can adapt the following code:

Code

costs.optimised <- c(
  NonSurgical = # Add the cost,
  SuccessfulNonSurgical = # Add the cost,
  TKR = # Add the cost,
  SuccessfulTKR = # Add the cost,
  Revision = # Add the cost,
  SuccessfulRevision = # Add the cost,
  Death = # Add the cost
)

Inspect the two vectors and confirm that their names appear in the same order as state.names:

Code

costs.current
costs.optimised

names(costs.current)
state.names

You can formally check the order using:

Code

identical(
  names(costs.current),
  state.names
)

identical(
  names(costs.optimised),
  state.names
)

Both results should be TRUE.

The order must match

Matrix multiplication uses the position of each value, not its meaning. If the order of the cost vector does not match the columns of the cohort trace, the code may still run but assign the wrong costs to the health states.

R will happily perform the calculation. Whether it is the calculation you intended is, regrettably, a separate question.

2. Create the state-utility vectors

Create a vector containing the utility values for current care:

Code

utilities.current <- c(
  NonSurgical = 0.658,
  SuccessfulNonSurgical = 0.749,
  TKR = 0.658,
  SuccessfulTKR = 0.878,
  Revision = 0.658,
  SuccessfulRevision = 0.760,
  Death = 0
)

Next, create a vector named utilities.optimised. The only difference is the utility of SuccessfulNonSurgical, which is 0.783 under optimised care.

Code

utilities.optimised <- c(
  NonSurgical = # Add the utility,
  SuccessfulNonSurgical = # Add the utility,
  TKR = # Add the utility,
  SuccessfulTKR = # Add the utility,
  Revision = # Add the utility,
  SuccessfulRevision = # Add the utility,
  Death = # Add the utility
)

Check that the names and order of both utility vectors match state.names:

Code

identical(
  names(utilities.current),
  state.names
)

identical(
  names(utilities.optimised),
  state.names
)

3. Calculate costs during each cycle

The cohort traces contain 41 rows, but the model contains only 40 cycles. The row for age 90 represents the distribution after the final transition and does not begin another cycle.

We therefore use rows 1 to n.cycles when calculating costs. Use:

Code

cycle.costs.current <- as.numeric(
  cohort.current[1:n.cycles, ] %*%
    costs.current
)

This calculation multiplies the number of people in each state by the corresponding state cost and sums the results across all states.

For example, the costs in cycle $t$ are:

\[ C_t = \sum_{j=1}^{7}m_{tj}c_j, \]

where:

$m_{tj}$ is the number of people in state $j$ during cycle $t$;
$c_j$ is the cost assigned to state $j$.

Adapt the code above to also calculate the cycle costs for optimised care and store them as:

cycle.costs.optimised.

4. Calculate QALYs during each cycle

Calculate QALYs under current care:

Code

cycle.qalys.current <- as.numeric(
  cohort.current[1:n.cycles, ] %*%
    utilities.current
)

Because one cycle represents one year, multiplying the number of people in each state by its utility gives the QALYs generated during that cycle.

Adapt the code to calculate QALYs under optimised care and store them as:

cycle.qalys.optimised.

5. Combine the cycle-specific results

Create a data frame containing the undiscounted cycle-specific results for current care:

Code

outcomes.current <- data.frame(
  cycle = 1:(n.cycles),
  age = (start.age + 1):(start.age + n.cycles),
  costs = cycle.costs.current,
  qalys = cycle.qalys.current
)

Create a corresponding data frame named outcomes.optimised using the optimised-care costs and QALYs.

Note: Cycle 0 is the initial state of the cohort before any model cycle has been completed. It is required in the cohort trace, but it should not be used to label the cycle-specific outcome results.

Inspect the first six rows:

Code

head(outcomes.current)

head(outcomes.optimised)

Compare the early cycles and consider:

why optimised care is more expensive in the first cycle;
why the QALYs differ after participants enter SuccessfulNonSurgical;
how delayed knee replacement affects costs in later cycles;
why the cycle-specific results represent totals for the complete cohort rather than values per person.

Do not sum the costs or QALYs yet. They are currently undiscounted, and the values occurring in later cycles must first be adjusted for time preference.

1. Create the state-cost vectors

We first create the cost vector for current care:

Code

costs.current <- c(
  NonSurgical = 222,
  SuccessfulNonSurgical = 200,
  TKR = 18326,
  SuccessfulTKR = 500,
  Revision = 28776,
  SuccessfulRevision = 500,
  Death = 0
)

costs.current

          NonSurgical SuccessfulNonSurgical                   TKR 
                  222                   200                 18326 
        SuccessfulTKR              Revision    SuccessfulRevision 
                  500                 28776                   500 
                Death 
                    0

We then create the corresponding vector for optimised care:

Code

costs.optimised <- c(
  NonSurgical = 1209,
  SuccessfulNonSurgical = 200,
  TKR = 18326,
  SuccessfulTKR = 500,
  Revision = 28776,
  SuccessfulRevision = 500,
  Death = 0
)

costs.optimised

          NonSurgical SuccessfulNonSurgical                   TKR 
                 1209                   200                 18326 
        SuccessfulTKR              Revision    SuccessfulRevision 
                  500                 28776                   500 
                Death 
                    0

The strategies have different costs only in the initial NonSurgical state. Optimised non-surgical treatment costs CHF 1,209 per person, compared with CHF 222 under current care. All subsequent state costs are identical.

We check whether the names and order match state.names:

Code

identical(
  names(costs.current),
  state.names
)

[1] TRUE

Code

identical(
  names(costs.optimised),
  state.names
)

[1] TRUE

Both results are TRUE. The state costs will therefore be multiplied by the correct columns of the cohort traces.

2. Create the state-utility vectors

We create the utility vector for current care:

Code

utilities.current <- c(
  NonSurgical = 0.658,
  SuccessfulNonSurgical = 0.749,
  TKR = 0.658,
  SuccessfulTKR = 0.878,
  Revision = 0.658,
  SuccessfulRevision = 0.760,
  Death = 0
)

utilities.current

          NonSurgical SuccessfulNonSurgical                   TKR 
                0.658                 0.749                 0.658 
        SuccessfulTKR              Revision    SuccessfulRevision 
                0.878                 0.658                 0.760 
                Death 
                0.000

We then create the utility vector for optimised care:

Code

utilities.optimised <- c(
  NonSurgical = 0.658,
  SuccessfulNonSurgical = 0.783,
  TKR = 0.658,
  SuccessfulTKR = 0.878,
  Revision = 0.658,
  SuccessfulRevision = 0.760,
  Death = 0
)

utilities.optimised

          NonSurgical SuccessfulNonSurgical                   TKR 
                0.658                 0.783                 0.658 
        SuccessfulTKR              Revision    SuccessfulRevision 
                0.878                 0.658                 0.760 
                Death 
                0.000

The strategies differ only in the utility assigned to SuccessfulNonSurgical. The higher value under optimised care represents the larger improvement in health-related quality of life produced by guideline-recommended non-surgical treatment.

We again check the order of the utility vectors:

Code

identical(
  names(utilities.current),
  state.names
)

[1] TRUE

Code

identical(
  names(utilities.optimised),
  state.names
)

[1] TRUE

Both results are again TRUE.

3. Calculate costs during each cycle

We calculate the costs generated during each cycle under current care:

Code

cycle.costs.current <- as.numeric(
  cohort.current[1:n.cycles, ] %*%
    costs.current
)

cycle.costs.current

 [1]  88800.00  79882.64 137664.31 137882.27 139473.61 141016.06 259612.33
 [8] 258394.62 259229.73 259947.22 260464.54 260824.63 260962.90 260878.07
[15] 260530.20 259916.64 374576.69 364895.95 357146.33 349255.90 341339.21
[22] 333522.32 325344.28 317104.76 308629.60 299872.68 300038.84 288020.86
[29] 276060.53 263885.18 251145.35 238376.01 225103.11 211187.99 197001.00
[36] 182384.04 123879.39 112696.41 100789.09  88552.33

The expression:

Code

cohort.current[1:n.cycles, ] %*%
  costs.current

multiplies the number of people in each health state by the corresponding state cost and sums the results across all states.

The function as.numeric() converts the resulting one-column matrix into a numeric vector.

We use only the first 40 rows of the cohort trace. Row 41 contains the cohort distribution at age 90 after the final transition and does not begin another model cycle.

We calculate the corresponding costs under optimised care:

Code

cycle.costs.optimised <- as.numeric(
  cohort.optimised[1:n.cycles, ] %*%
    costs.optimised
)

cycle.costs.optimised

 [1] 483600.00  79882.64  79746.20  79569.96 137084.64 137267.43 258556.72
 [8] 257358.86 258247.64 259018.09 259586.20 259996.17 260182.69 260145.19
[15] 259843.11 259274.48 376381.19 366587.19 358758.23 350791.51 342798.04
[22] 334912.47 326663.68 318358.11 309816.91 300998.72 301589.78 289459.96
[29] 277390.79 265118.04 252277.51 239414.91 226052.93 212050.84 197778.47
[36] 183082.86 123956.77 112773.08 100861.18  88619.08

In the first cycle, all 400 people occupy NonSurgical. The current-care costs are therefore:

\[ 400\times CHF\ 222 = CHF\ 88{,}800. \]

Under optimised care, the corresponding costs are:

\[ 400\times CHF\ 1{,}209 = CHF\ 483{,}600. \]

The optimised strategy is therefore considerably more expensive in the first cycle. In later cycles, differences arise from the number of people occupying the surgical and post-surgical states.

4. Calculate QALYs during each cycle

We calculate QALYs under current care:

Code

cycle.qalys.current <- as.numeric(
  cohort.current[1:n.cycles, ] %*%
    utilities.current
)

cycle.qalys.current

 [1] 263.2000 299.1605 298.3587 298.1128 297.8494 297.5621 296.5922 296.8694
 [9] 297.1277 297.1584 297.1090 296.8354 296.3944 295.6787 294.7478 293.4930
[17] 291.4373 290.8868 289.9126 288.3839 286.8887 284.5529 281.9877 278.7641
[25] 275.2360 270.7173 265.9957 260.4046 254.7196 247.3836 239.6989 231.1366
[33] 221.3737 210.4667 199.0705 186.2371 171.3498 155.2119 138.3289 121.0420

We then calculate QALYs under optimised care:

Code

cycle.qalys.optimised <- as.numeric(
  cohort.optimised[1:n.cycles, ] %*%
    utilities.optimised
)

cycle.qalys.optimised

 [1] 263.2000 312.7405 312.2064 311.5164 310.4267 310.0268 308.6962 308.6407
 [9] 308.5768 308.2874 307.9255 307.3419 306.5953 305.5750 304.3426 302.7863
[17] 300.1831 299.1429 297.6980 295.7137 293.7928 291.0383 288.0764 284.4676
[25] 280.5735 275.6938 270.5877 264.6274 258.6009 250.9254 242.9240 234.0594
[33] 224.0043 212.8166 201.1590 188.0723 173.0010 156.6720 139.5987 122.1259

In the first cycle, all 400 people occupy NonSurgical, which has a utility of 0.658 under both strategies. The cohort therefore generates:

\[ 400\times0.658 = 263.2\text{ QALYs}. \]

The first-cycle QALYs are identical because the strategies assign the same utility to NonSurgical.

From the second cycle onwards, most surviving participants occupy SuccessfulNonSurgical. The utility of this state is higher under optimised care:

\[ 0.783>0.749. \]

Optimised care therefore generates more QALYs while participants remain in successful non-surgical treatment.

5. Combine the cycle-specific results

We combine the current-care results in one data frame:

Code

outcomes.current <- data.frame(
  cycle = 1:(n.cycles),
  age = (start.age + 1):(start.age + n.cycles),
  costs = cycle.costs.current,
  qalys = cycle.qalys.current
)

We create the corresponding data frame for optimised care:

Code

outcomes.optimised <- data.frame(
  cycle = 1:(n.cycles),
  age = (start.age + 1):(start.age + n.cycles),
  costs = cycle.costs.optimised,
  qalys = cycle.qalys.optimised
)

We inspect the first six cycles:

Code

head(outcomes.current)

  cycle age     costs    qalys
1     1  51  88800.00 263.2000
2     2  52  79882.64 299.1605
3     3  53 137664.31 298.3587
4     4  54 137882.27 298.1128
5     5  55 139473.61 297.8494
6     6  56 141016.06 297.5621

Code

head(outcomes.optimised)

  cycle age     costs    qalys
1     1  51 483600.00 263.2000
2     2  52  79882.64 312.7405
3     3  53  79746.20 312.2064
4     4  54  79569.96 311.5164
5     5  55 137084.64 310.4267
6     6  56 137267.43 310.0268

Optimised care is more expensive in cycle 0 because all participants receive the more intensive initial treatment.

In cycle 1, costs are the same under both strategies because all surviving participants occupy SuccessfulNonSurgical, which costs CHF 200 under both strategies. However, optimised care generates more QALYs because it assigns a higher utility to this state.

Under current care, total knee replacements first generate costs in cycle 2. Under optimised care, these costs are delayed because no participants enter TKR during the transitions beginning at ages 51 and 52.

The cycle-specific values are totals for the complete cohort. They have not yet been discounted or converted into totals per person.

3.4 4. Calculate total discounted costs and QALYs

The previous step calculated the costs and QALYs generated during each model year/cycle. The first model year extends from age 50 to 51, while the final model year extends from age 89 to 90. For discounting, these years are indexed from $t=0$ to $t=39$. In the context of discounting, the first model year is assigned $t=0$ because its costs and QALYs are treated as occurring without delay and are therefore not discounted.

Health economic evaluations generally give less weight to costs and health outcomes that occur further in the future. This reflects time preference: costs and health benefits occurring today are usually valued more highly than equivalent costs and benefits occurring many years later.

We therefore convert future costs and QALYs into their present values using an annual discount rate of 3%.

For a value occurring in year $t$, the discount factor is:

\[ DF_t = \frac{1}{(1+r)^t}, \]

where:

$DF_t$ is the discount factor for year $t$;
$r$ is the annual discount rate;
$t$ is the number of years after the beginning of the model.

The discounted cost in year $t$ is:

\[ PV(C_t) = C_t\times DF_t = \frac{C_t}{(1+r)^t}, \]

where:

$C_t$ is the undiscounted cost generated during year $t$;
$PV(C_t)$ is its present value.

The discounted QALYs are calculated in the same way:

\[ PV(Q_t) = Q_t\times DF_t = \frac{Q_t}{(1+r)^t}, \]

where:

$Q_t$ is the undiscounted number of QALYs generated during year $t$;
$PV(Q_t)$ is the present value of those QALYs.

As the first year belongs to $t = 0$, it is not discounted:

\[ DF_0 = \frac{1}{(1.03)^0} = 1. \]

The values generated during year 2 ($t = 1$) are discounted by one year:

\[ DF_1 = \frac{1}{(1.03)^1} \approx 0.9709. \]

For example, a cost of CHF 10,000 occurring in year 2 has a present value of:

\[ PV(C_1) = \frac{CHF\ 10{,}000}{1.03} \approx CHF\ 9{,}709. \]

A cost occurring in year 21 ($t = 20$) receives substantially less weight:

\[ DF_{20} = \frac{1}{(1.03)^{20}} \approx 0.5537. \]

The discounted total costs over all 40 model years are:

\[ C_{\text{total}} = \sum_{t=0}^{39} \frac{C_t}{(1+r)^t}, \]

and the discounted total QALYs are:

\[ Q_{\text{total}} = \sum_{t=0}^{39} \frac{Q_t}{(1+r)^t}. \]

These calculations initially produce totals for the complete cohort of 400 people. To obtain the expected costs and QALYs per person, we divide the cohort totals by the initial cohort size:

\[ C_{\text{per person}} = \frac{C_{\text{total}}}{400}, \]

and:

\[ Q_{\text{per person}} = \frac{Q_{\text{total}}}{400}. \]

Per-person values allow the two strategies to be compared and will later be used to calculate incremental costs, incremental QALYs and the incremental cost-effectiveness ratio.

Discounting is not cost standardisation

The model costs are expressed in 2019 Swiss francs. Discounting does not convert them into another price year.

Cost standardisation adjusts monetary values for changes in prices between calendar years.
Discounting adjusts future costs and outcomes to reflect when they occur.

The calculations may look similar because both involve adjustment factors, but they answer different questions. Health economics enjoys keeping these distinctions clear, even when the formulas attempt to look related.

Exercise

1. Create the discount factors

Set the annual discount rate to 3%:

Code

discount.rate <- 0.03

Calculate one discount factor for every model year/cycle. Use:

Code

discount.factor <- 1 /
  (1 + discount.rate) ^
  (outcomes.current$cycle - 1)

Inspect the first and final discount factors:

Code

head(discount.factor)

tail(discount.factor)

Check that:

the discount factor for cycle 0 is 1;
the discount factors become smaller over time;
the vector contains 40 values.

You can check its length using:

Code

length(discount.factor)

2. Discount the cycle-specific costs

Add the discount factors to the current-care outcomes:

Code

outcomes.current$discount.factor <-
  discount.factor

Calculate the present value of the costs generated during every cycle:

Code

outcomes.current$discounted.costs <-
  outcomes.current$costs *
  outcomes.current$discount.factor

Add the same discount factors to outcomes.optimised and calculate:

outcomes.optimised$discounted.costs.

Both strategies must use the same discount factors because their costs occur over the same model cycles.

Inspect the first rows of the current-care results:

Code

head(outcomes.current)

Confirm that costs in cycle 0 remain unchanged and that costs in later cycles are reduced.

3. Discount the cycle-specific QALYs

Calculate the discounted QALYs under current care:

Code

outcomes.current$discounted.qalys <-
  outcomes.current$qalys *
  outcomes.current$discount.factor

Adapt this code to calculate:

outcomes.optimised$discounted.qalys.

Inspect the undiscounted and discounted QALYs:

Code

head(outcomes.current)

As with costs, QALYs generated during cycle 0 remain unchanged, while QALYs generated in later cycles receive progressively less weight.

4. Calculate total discounted costs and QALYs

Calculate the total discounted costs for the complete current-care cohort:

Code

total.costs.current <- sum(
  outcomes.current$discounted.costs
)

Calculate the total discounted QALYs:

Code

total.qalys.current <- sum(
  outcomes.current$discounted.qalys
)

Adapt these calculations for optimised care and store the results as:

total.costs.optimised;
total.qalys.optimised.

These values are totals for the complete cohort of 400 people.

5. Calculate costs and QALYs per person

Calculate the initial cohort size rather than entering the number 400 again:

Code

cohort.size <- sum(
  cohort.initial
)

Calculate discounted costs and QALYs per person under current care:

Code

costs.per.person.current <-
  total.costs.current /
  cohort.size

qalys.per.person.current <-
  total.qalys.current /
  cohort.size

Adapt these calculations for optimised care and store the results as:

costs.per.person.optimised;
qalys.per.person.optimised.

Display the four per-person results:

Code

costs.per.person.current
costs.per.person.optimised

qalys.per.person.current
qalys.per.person.optimised

Do not yet calculate incremental costs, incremental QALYs or the ICER. These comparisons follow in the next step.

6. Check the effect of discounting

Compare the discounted and undiscounted totals for current care:

Code

sum(outcomes.current$costs)
total.costs.current

sum(outcomes.current$qalys)
total.qalys.current

Repeat the comparison for optimised care.

Confirm that:

discounted totals are lower than the corresponding undiscounted totals;
cycle 1 is not discounted;
the same discount factors are used for both strategies;
the per-person values equal the cohort totals divided by the initial cohort size.

Keep the levels of aggregation clear

The columns in outcomes.current and outcomes.optimised contain results for the complete cohort during each cycle.

The objects beginning with total. contain results for the complete cohort over all 40 cycles.

The objects beginning with costs.per.person or qalys.per.person contain the expected results for one person over the complete modelled period.

Keeping these levels separate prevents a very easy mistake: comparing the costs of 400 people with the QALYs of one person. That would certainly produce a ratio, although not one we would wish to defend.

1. Create the discount factors

We first define the annual discount rate:

Code

discount.rate <- 0.03

We then calculate one discount factor for each of the 40 model cycles:

Code

discount.factor <- 1 /
  (1 + discount.rate) ^
  (outcomes.current$cycle - 1)

We inspect the first values:

Code

head(discount.factor)

[1] 1.0000000 0.9708738 0.9425959 0.9151417 0.8884870 0.8626088

Cycle 0 has a discount factor of 1 and therefore remains undiscounted. The factors become progressively smaller because values occurring further in the future receive less weight.

We inspect the final values:

Code

tail(discount.factor)

[1] 0.3660449 0.3553834 0.3450324 0.3349829 0.3252262 0.3157535

We also check the length of the vector:

Code

length(discount.factor)

[1] 40

The result is 40, corresponding to the 40 cycles from age 50–51 to age 89–90.

2. Discount the cycle-specific costs

We add the discount factors to the current-care outcomes:

Code

outcomes.current$discount.factor <-
  discount.factor

We calculate the present value of the costs generated during each cycle:

Code

outcomes.current$discounted.costs <-
  outcomes.current$costs *
  outcomes.current$discount.factor

We perform the same calculations for optimised care:

Code

outcomes.optimised$discount.factor <-
  discount.factor

outcomes.optimised$discounted.costs <-
  outcomes.optimised$costs *
  outcomes.optimised$discount.factor

Both strategies use the same discount factors because they cover the same model cycles.

We inspect the first current-care values:

Code

head(outcomes.current)

  cycle age     costs    qalys discount.factor discounted.costs
1     1  51  88800.00 263.2000       1.0000000         88800.00
2     2  52  79882.64 299.1605       0.9708738         77555.96
3     3  53 137664.31 298.3587       0.9425959        129761.82
4     4  54 137882.27 298.1128       0.9151417        126181.81
5     5  55 139473.61 297.8494       0.8884870        123920.49
6     6  56 141016.06 297.5621       0.8626088        121641.69

The first-cycle costs therefore remain unchanged. From cycle 2 onwards, discounted costs are lower than the corresponding undiscounted costs.

3. Discount the cycle-specific QALYs

We calculate discounted QALYs under current care:

Code

outcomes.current$discounted.qalys <-
  outcomes.current$qalys *
  outcomes.current$discount.factor

We repeat the calculation for optimised care:

Code

outcomes.optimised$discounted.qalys <-
  outcomes.optimised$qalys *
  outcomes.optimised$discount.factor

We inspect the first current-care values:

Code

head(outcomes.current)

  cycle age     costs    qalys discount.factor discounted.costs
1     1  51  88800.00 263.2000       1.0000000         88800.00
2     2  52  79882.64 299.1605       0.9708738         77555.96
3     3  53 137664.31 298.3587       0.9425959        129761.82
4     4  54 137882.27 298.1128       0.9151417        126181.81
5     5  55 139473.61 297.8494       0.8884870        123920.49
6     6  56 141016.06 297.5621       0.8626088        121641.69
  discounted.qalys
1         263.2000
2         290.4471
3         281.2317
4         272.8155
5         264.6353
6         256.6797

QALYs generated during cycle 1 remain unchanged. QALYs generated during later cycles receive progressively less weight.

4. Calculate total discounted costs and QALYs

We sum the discounted costs under current care:

Code

total.costs.current <- sum(
  outcomes.current$discounted.costs
)

total.costs.current

[1] 5583624

We sum the discounted QALYs:

Code

total.qalys.current <- sum(
  outcomes.current$discounted.qalys
)

total.qalys.current

[1] 6513.129

We repeat the calculations for optimised care:

Code

total.costs.optimised <- sum(
  outcomes.optimised$discounted.costs
)

total.qalys.optimised <- sum(
  outcomes.optimised$discounted.qalys
)

We display the results:

Code

round(
  c(
    total.costs.current = total.costs.current,
    total.costs.optimised = total.costs.optimised,
    total.qalys.current = total.qalys.current,
    total.qalys.optimised = total.qalys.optimised
  ),
  2
)

  total.costs.current total.costs.optimised   total.qalys.current 
           5583623.91            5871380.39               6513.13 
total.qalys.optimised 
              6705.12

These values represent the complete cohort of 400 people over the 40-year modelled period.

5. Calculate costs and QALYs per person

We calculate the initial cohort size:

Code

cohort.size <- sum(
  cohort.initial
)

cohort.size

[1] 400

The result is 400.

We calculate costs and QALYs per person under current care:

Code

costs.per.person.current <-
  total.costs.current /
  cohort.size

qalys.per.person.current <-
  total.qalys.current /
  cohort.size

We repeat the calculations for optimised care:

Code

costs.per.person.optimised <-
  total.costs.optimised /
  cohort.size

qalys.per.person.optimised <-
  total.qalys.optimised /
  cohort.size

We display the per-person results:

Code

round(
  c(
    costs.per.person.current =
      costs.per.person.current,
    costs.per.person.optimised =
      costs.per.person.optimised,
    qalys.per.person.current =
      qalys.per.person.current,
    qalys.per.person.optimised =
      qalys.per.person.optimised
  ),
  2
)

  costs.per.person.current costs.per.person.optimised 
                  13959.06                   14678.45 
  qalys.per.person.current qalys.per.person.optimised 
                     16.28                      16.76

6. Check the effect of discounting

We compare discounted and undiscounted costs under current care:

Code

sum(outcomes.current$costs)

[1] 9556288

Code

total.costs.current

[1] 5583624

We compare the corresponding QALYs:

Code

sum(outcomes.current$qalys)

[1] 10457.44

Code

total.qalys.current

[1] 6513.129

We repeat the comparison for optimised care:

Code

sum(outcomes.optimised$costs)

[1] 9845852

Code

total.costs.optimised

[1] 5871380

Code

sum(outcomes.optimised$qalys)

[1] 10734.44

Code

total.qalys.optimised

[1] 6705.118

All discounted totals are lower than the corresponding undiscounted totals because future values receive less weight.

The difference becomes substantial over a 40-year model because costs and QALYs arising near the end of the time horizon are discounted for many years.

3.5 5. Compare current and optimised care

The final step is to compare the expected costs and QALYs of the two strategies. Health economic evaluation focuses on the incremental results: the additional costs and additional health outcomes produced by optimised care relative to current care.

The incremental cost per person is:

\[ \Delta C = C_{\text{optimised}} - C_{\text{current}}, \]

where:

$\Delta C$ is the incremental cost;
$C_{\text{optimised}}$ is the discounted cost per person under optimised care;
$C_{\text{current}}$ is the discounted cost per person under current care.

The incremental QALYs are:

\[ \Delta Q = Q_{\text{optimised}} - Q_{\text{current}}, \]

where:

$\Delta Q$ is the incremental number of QALYs;
$Q_{\text{optimised}}$ is the discounted number of QALYs per person under optimised care;
$Q_{\text{current}}$ is the discounted number of QALYs per person under current care.

The signs of these differences should be examined before calculating a ratio:

Incremental costs	Incremental QALYs	Interpretation
Positive	Positive	Optimised care is more effective and more costly
Negative	Positive	Optimised care is dominant
Positive	Negative	Optimised care is dominated
Negative	Negative	Optimised care is less effective and less costly

If optimised care is more effective and more costly, we calculate the incremental cost-utility ratio:

\[ ICUR = \frac{\Delta C}{\Delta Q}. \]

The ICUR represents the additional cost required to gain one additional QALY with optimised care.

For example, an ICUR of CHF 20,000 per QALY means that optimised care generates additional health benefits at an additional cost of CHF 20,000 for each additional QALY gained.

Whether this represents good value depends on the decision-maker’s willingness to pay for an additional QALY. Following the published evaluation, we use a willingness-to-pay threshold of CHF 100,000 per QALY.

If:

\[ ICUR<CHF\ 100{,}000\text{ per QALY}, \]

optimised care is considered cost-effective under this decision rule.

The same decision can also be expressed using incremental net monetary benefit:

\[ INMB = \lambda\Delta Q-\Delta C, \]

where:

$INMB$ is the incremental net monetary benefit;
$\lambda$ is the willingness-to-pay threshold;
$\Delta Q$ is the incremental QALY gain;
$\Delta C$ is the incremental cost.

If:

\[ INMB>0, \]

the monetary value assigned to the additional QALYs exceeds the additional cost, and optimised care is considered cost-effective at the selected threshold.

Net monetary benefit is especially helpful when incremental cost-effectiveness ratios are negative or difficult to interpret. Ratios may become rather dramatic when the denominator approaches zero; subtraction generally behaves itself rather better.

Interpretation applies to this model

The result represents the expected cost-effectiveness of optimised care under the structure, parameters and assumptions used in this exercise.

It does not demonstrate that optimised care will delay every knee replacement by exactly two years in clinical practice. In particular, the result depends on:

the assumed two-year delay;
the state costs and utilities;
the age-specific transition probabilities;
the 40-year time horizon;
the 3% discount rate;
the simplified modelling of mortality and transitions.

A favourable ICUR is therefore a model-based estimate, not a promise from the future.

Exercise

1. Summarise the per-person results

Create a data frame containing the discounted costs and QALYs per person under both strategies:

Code

strategy.results <- data.frame(
  strategy = c(
    "Current care",
    "Optimised care"
  ),
  costs = c(
    costs.per.person.current,
    costs.per.person.optimised
  ),
  qalys = c(
    qalys.per.person.current,
    qalys.per.person.optimised
  )
)

strategy.results

Have a look at strategy.results.

2. Calculate incremental costs and QALYs

Calculate the incremental cost of optimised care:

Code

delta.cost <-
  costs.per.person.optimised -
  costs.per.person.current

Calculate the incremental QALYs and store the result as delta.qaly.

Display both results:

Code

delta.cost
delta.qaly

Interpret their signs:

Is optimised care more or less costly than current care?
Does optimised care generate more or fewer QALYs?
In which quadrant of the cost-effectiveness plane does the result belong?

3. Calculate the ICUR

Calculate the incremental cost-utility ratio:

Code

icur <-
  delta.cost /
  delta.qaly

icur

The ICUR represents the additional cost per additional QALY gained with optimised care.

Round the result to the nearest Swiss franc:

Code

round(icur, 0)

Before interpreting the ICUR, check that both delta.cost and delta.qaly are positive. If either value is negative, the ratio requires a different interpretation.

4. Compare the ICUR with the willingness-to-pay threshold

Set the willingness-to-pay threshold to CHF 100,000 per QALY:

Code

wtp <- 100000

Use a logical comparison to determine whether the ICUR is below the threshold:

Code

icur < wtp

A result of TRUE indicates that optimised care is considered cost-effective according to this decision rule.

5. Calculate incremental net monetary benefit

Calculate the incremental net monetary benefit at the same willingness-to-pay threshold:

Code

inmb <-
  wtp * delta.qaly -
  delta.cost

inmb

Interpret the result:

inmb > 0 supports optimised care;
inmb < 0 supports current care;
inmb = 0 means that the strategies have the same net monetary benefit at this threshold.

Confirm the decision using:

Code

inmb > 0

The conclusions based on the ICUR and incremental net monetary benefit should agree.

6. Report the result

Write a short interpretation that includes:

the discounted cost per person under each strategy;
the discounted QALYs per person under each strategy;
the incremental cost;
the incremental QALY gain;
the ICUR;
the willingness-to-pay threshold;
the incremental net monetary benefit;
the resulting decision.

Cost-effective does not mean cost-saving

An intervention is cost-saving or dominant, when it produces better outcomes at lower costs.

An intervention is cost-effective when its additional health benefits are considered sufficient to justify its additional costs at the selected willingness-to-pay threshold.

Optimised care can therefore be cost-effective without reducing total healthcare costs. Better health sometimes costs more; the economic question is whether the additional benefit is worth the additional cost.

1. Summarise the per-person results

We combine the discounted per-person results in one data frame:

Code

strategy.results <- data.frame(
  strategy = c(
    "Current care",
    "Optimised care"
  ),
  costs = c(
    costs.per.person.current,
    costs.per.person.optimised
  ),
  qalys = c(
    qalys.per.person.current,
    qalys.per.person.optimised
  )
)

strategy.results

        strategy    costs    qalys
1   Current care 13959.06 16.28282
2 Optimised care 14678.45 16.76280

Optimised care produces more QALYs but also has higher expected costs.

2. Calculate incremental costs and QALYs

We calculate the incremental cost:

Code

delta.cost <-
  costs.per.person.optimised -
  costs.per.person.current

delta.cost

[1] 719.3912

Optimised care therefore costs approximately CHF 719 more per person than current care.

We calculate the incremental QALYs:

Code

delta.qaly <-
  qalys.per.person.optimised -
  qalys.per.person.current

delta.qaly

[1] 0.4799743

Optimised care therefore generates approximately 0.48 additional discounted QALYs per person.

Both incremental costs and incremental QALYs are positive. The result lies in the north-east quadrant of the cost-effectiveness plane:

optimised care is more effective;
optimised care is more costly.

Whether it is cost-effective therefore depends on the willingness to pay for the additional QALYs.

3. Calculate the ICUR

We calculate the incremental cost-utility ratio:

Code

icur <-
  delta.cost /
  delta.qaly

icur

[1] 1498.812

Rounded to the nearest Swiss franc:

Code

round(icur, 0)

[1] 1499

The ICUR is approximately CHF 1,499 per additional QALY gained.

Both the numerator and denominator are positive, so the ratio has a straightforward interpretation: optimised care generates additional QALYs at an additional cost of approximately CHF 1,499 per QALY.

4. Compare the ICUR with the willingness-to-pay threshold

We define the willingness-to-pay threshold:

Code

wtp <- 100000

We compare the ICUR with the threshold:

Code

icur < wtp

[1] TRUE

The ICUR of approximately CHF 1,499 per QALY is substantially below the willingness-to-pay threshold of CHF 100,000 per QALY.

According to this decision rule, optimised care is cost-effective compared with current care.

5. Calculate incremental net monetary benefit

We calculate the incremental net monetary benefit:

Code

inmb <-
  wtp * delta.qaly -
  delta.cost

inmb

[1] 47278.04

At a willingness-to-pay threshold of CHF 100,000 per QALY, the monetary value assigned to the additional QALYs exceeds the additional costs by approximately CHF 47,278 per person.

The incremental net monetary benefit therefore supports optimised care. This conclusion agrees with the interpretation of the ICUR.

6. Report the result

Over the 40-year modelled period, current care produced approximately 16.28 discounted QALYs at a discounted cost of CHF 13,959 per person. Optimised care produced approximately 16.76 discounted QALYs at a cost of CHF 14,678 per person.

Compared with current care, optimised care generated approximately 0.48 additional QALYs at an additional cost of CHF 719 per person. This resulted in an ICUR of approximately CHF 1,499 per additional QALY.

The ICUR was below the willingness-to-pay threshold of CHF 100,000 per QALY. At this threshold, the incremental net monetary benefit was approximately CHF 47,278 per person. Under the assumptions of this model, optimised care is therefore considered cost-effective but not cost-saving.

Simplifications of the teaching model

The results differ from those reported in the original publication because this teaching model uses a simplified structure and modified assumptions. The purpose is to demonstrate the main steps of Markov modelling, not to reproduce the published results exactly.

3.6 Sensitivity analysis

The results of the Markov model depend on uncertain input parameters and modelling assumptions. The base-case result should therefore not be interpreted as an exact prediction. Sensitivity analysis examines whether the conclusion changes when these inputs or assumptions are varied.

A probabilistic sensitivity analysis would be particularly appropriate because several parameters are uncertain simultaneously, including:

the costs of non-surgical treatment, TKR and revision surgery;
the utility values assigned to the health states;
the probabilities of TKR and revision surgery.

In a probabilistic sensitivity analysis, uncertain parameters are assigned probability distributions. For example:

beta distributions can be used for probabilities and utilities because their values are bounded;
gamma distributions can be used for costs because costs are positive and often right-skewed.

One value is randomly drawn from each distribution, the complete Markov model is recalculated, and the incremental costs and QALYs are stored. This procedure is repeated many times:

\[ \left( \Delta C^{(b)}, \Delta Q^{(b)} \right), \qquad b=1,\ldots,B, \]

where:

$b$ identifies one probabilistic simulation;
$B$ is the total number of simulations;
$\Delta C^{(b)}$ is the incremental cost in simulation $b$;
$\Delta Q^{(b)}$ is the incremental QALY gain in simulation $b$.

The simulated incremental costs and QALYs can be displayed on a cost-effectiveness plane. They can also be used to estimate the probability that optimised care is cost-effective at different willingness-to-pay thresholds:

\[ \Pr\left( \lambda\Delta Q-\Delta C>0 \right). \]

The assumed two-year delay of TKR is slightly different. It is a structural assumption, not simply an uncertain numerical parameter. It would therefore be more appropriate to examine it using scenario analysis, for example by comparing:

no delay;
a two-year delay;
a five-year delay.

A complete probabilistic sensitivity analysis would require probability distributions and measures of uncertainty for all relevant parameters. Because these inputs and the additional programming would considerably expand the exercise, we do not implement it here. Nevertheless, it would be an important part of a publication-quality model.

4 References and further reading

The following references provide methodological guidance and applied examples relevant to the trial-based and model-based health economic evaluations covered in this exercise.

4.1 General guidance and reporting

Husereau D, Drummond M, Augustovski F, et al. Consolidated Health Economic Evaluation Reporting Standards 2022 (CHEERS 2022) statement: updated reporting guidance for health economic evaluations. BMJ. 2022;376:e067975.
The CHEERS checklist provides guidance on the transparent reporting of trial-based and model-based economic evaluations.
Annemans L. Health economics for non-economists: principles, methods and pitfalls of health economic evaluations. 2nd ed. Kalmthout: Pelckmans Pro; 2018.
Taeymans J, Pfeiffer F. Gesundheitsökonomische Evaluationen physiotherapeutischer Interventionen. physioscience. 2017;13(1):9–16.

4.2 Trial-based economic evaluation

Jornada Ben Â, van Dongen JM, El Alili M, Esser JL, Broulíková HM, Bosmans JE. Conducting trial-based economic evaluations using R: a tutorial. PharmacoEconomics. 2023;41(11):1403–1413.
This practical tutorial demonstrates how trial-based economic evaluations can be conducted in R. It discusses missing data, baseline imbalances, skewed costs and the correlation between costs and outcomes.
Ramsey SD, Willke RJ, Glick H, Reed SD, Augustovski F, Jonsson B, Briggs A, Sullivan SD. Cost-effectiveness analysis alongside clinical trials II: an ISPOR Good Research Practices Task Force report. Value in Health. 2015;18(2):161–172. This report provides methodological recommendations for planning, analysing and reporting economic evaluations conducted alongside clinical trials.
Manca A, Hawkins N, Sculpher MJ. Estimating mean QALYs in trial-based cost-effectiveness analysis: the importance of controlling for baseline utility. Health Economics. 2005;14(5):487–496.
This paper explains why baseline utility should be considered when comparing QALYs between treatment groups.

4.3 Markov and state-transition models

Sonnenberg FA, Beck JR. Markov models in medical decision making: a practical guide. Medical Decision Making. 1993;13(4):322–338.
This classic introduction explains health states, transition probabilities, cohort traces and the Markov assumption.
Roberts M, Russell LB, Paltiel AD, Chambers M, McEwan P, Krahn M. Conceptualizing a model: a report of the ISPOR-SMDM Modeling Good Research Practices Task Force–2. Medical Decision Making. 2012;32(5):678–689.
This paper discusses how a decision model should be designed so that its structure reflects the clinical problem and decision being examined.
Siebert U, Alagoz O, Bayoumi AM, Jahn B, Owens DK, Cohen DJ, Kuntz KM. State-transition modeling: a report of the ISPOR-SMDM Modeling Good Research Practices Task Force–3. Medical Decision Making. 2012;32(5):690–700.
This report provides good-practice guidance on constructing, analysing and validating state-transition models.
Briggs AH, Weinstein MC, Fenwick EAL, Karnon J, Sculpher MJ, Paltiel AD. Model parameter estimation and uncertainty: a report of the ISPOR-SMDM Modeling Good Research Practices Task Force–6. Value in Health. 2012;15(6):835–842.
This paper distinguishes parameter, structural and stochastic uncertainty and explains deterministic and probabilistic sensitivity analysis.
Carta A, Conversano C. On the use of Markov models in pharmacoeconomics: pros and cons and implications for policy makers. Frontiers in Public Health. 2020;8:569500.
This accessible overview discusses the strengths and limitations of Markov models in economic evaluation.

4.4 Applied examples

Schurz A, Deliens T, Liechti M, Vanroose M, Clijsen R, Nijs J, Malfliet A, Van Bogaert W, Clarys P, Baur H, Taeymans J, Lutz N. Economic evaluation of a lifestyle intervention for individuals with overweight or obesity suffering from chronic low back pain—the BO2WL trial: a protocol for a health economic analysis. BMJ Open. 2025;15(6):e098272.
This protocol illustrates how a trial-based economic evaluation can be planned alongside a physiotherapy-related intervention.
Lutz N, Deliens T, Clarys P, Verhaeghe N, Taeymans J. Health economic evaluation of an influenza vaccination programme to prevent sick leave in employees: a prospective cohort study. Journal of Occupational and Environmental Medicine. 2020;62(8):549–556.
This study provides an applied example in which healthcare costs and productivity outcomes are considered.
Vetsch T, Taeymans J, Lutz N. Optimising the current model of care for knee osteoarthritis with the implementation of guideline-recommended non-surgical treatments: a model-based health economic evaluation. Swiss Medical Weekly. 2023;153:40059.
This is the original Swiss Markov-model evaluation on which the model-based part of this exercise is based.