qalytools provides a simple and intuitive user interface for the analysis of EQ-5D surveys. It builds upon the eq5d package to facilitate the calculation of QALY metrics and other related values across multiple surveys.
Overview
The package provides a range of functions:
- Constructors (and validators) for EQ5D data frame subclasses (EQ5D3L, EQ5D5L and EQ5D):
new_eq5d3l()(validate_eq5d3l()),new_eq5d5l()(validate_eq5d5l()) andnew_eq5dy()(validate_eq5dy()). - User-friendly coercion functions
as_eq5d3l(),as_eq5d5l()andas_eq5dy(). - The calculation of utility values based on a range of different value sets. This functionality is provided via the
calculate_utility(),add_utility()andavailable_valuesets()functions which are wrappers around the eq5d package. - The calculation of different Quality of Life Years (QALY) metrics including unadjusted ‘raw’ values, and the disutility from both perfect health and, optionally, a specified baseline. See
calculate_qalys(). - The calculation of the Paretian Classification of Health Change (PCHC) in an individual’s health state between two surveys via
calculate_pchc()(again wrapping the eq5d package). - Easy calculation of responses with a health limitation (i.e. a non-one response in one of the dimensions) via
calculate_limitation().
Whilst this vignette provides an introduction to the core functionality of the package, a more detailed analysis example can be found in vignette("example_analysis").
The EQ5D object class
We define an EQ5D object as a table with data represented in long format that meets a few additional criteria:
It contains columns that represent dimensions from the EQ5D survey specification as well as a column representing the Visual Analogue Score (columns ‘mobility’, ‘self_care’, ‘usual’, ‘pain’, ‘anxiety’ and ‘vas’ in the example table below)
Dimension values must be whole numbers, bounded between 1 and 3 for EQ5D3L and EQ5D3Y surveys or bounded between 1 and 5 for EQ5D5L surveys
It contains a column that acts as a unique respondent identifier (
respondentID) and another that identifying different surveys over time (`surveyID). Together these should uniquely identify a response (i.e. no combination of these should be duplicated within the given data frame).
| surveyID | respondentID | time_index | mobility | self_care | usual | pain | anxiety | vas | sex | age |
|---|---|---|---|---|---|---|---|---|---|---|
| . | . | . | . | . | . | . | . | . | . | . |
| 1 | 1 | 30 | 3 | 2 | 2 | 1 | 1 | 87 | Female | 25 |
| 1 | 2 | 30 | 2 | 2 | 2 | 2 | 2 | 65 | Male | 30 |
| 1 | 3 | 30 | 1 | 3 | 3 | 1 | 3 | 76 | Male | 39 |
| 1 | 4 | 30 | 1 | 2 | 3 | 2 | 4 | 55 | Female | 60 |
| 1 | 5 | 30 | 2 | 1 | 4 | 1 | 1 | 80 | Male | 28 |
| 1 | 6 | 30 | 1 | 3 | 2 | 1 | 1 | 83 | Male | 59 |
| . | . | . | . | . | . | . | . | . | . | . |
In qalytools these EQ5D objects are implemented as a subclass of data frame and we provide associated methods for working with this class.
Usage
In the following examples we make use of a synthetic EQ-5D-5L data set that is included in the package
library(qalytools)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(ggplot2)
# Example EQ5D5L data
data("EQ5D5L_surveys")
glimpse(EQ5D5L_surveys)
#> Rows: 10,000
#> Columns: 12
#> $ surveyID <chr> "survey01", "survey02", "survey03", "survey04", "survey05…
#> $ respondentID <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, …
#> $ sex <chr> "Female", "Female", "Female", "Female", "Female", "Female…
#> $ age <int> 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 30, 30, 30, 30, 3…
#> $ mobility <int> 2, 1, 1, 1, 1, 2, 2, 1, 1, 2, 5, 5, 1, 1, 3, 1, 1, 2, 1, …
#> $ self_care <int> 1, 4, 5, 2, 2, 1, 1, 1, 3, 1, 1, 2, 3, 3, 3, 4, 4, 3, 2, …
#> $ usual <int> 2, 5, 2, 5, 3, 1, 2, 2, 2, 2, 1, 5, 5, 3, 5, 2, 2, 1, 3, …
#> $ pain <int> 1, 5, 1, 2, 3, 2, 1, 2, 1, 1, 1, 5, 1, 2, 1, 2, 2, 2, 3, …
#> $ anxiety <int> 1, 3, 3, 2, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 1, 1, 1, 2, 2, …
#> $ time_index <dbl> 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 30, 60, 90…
#> $ vas <dbl> 96, 22, 61, 65, 78, 88, 85, 95, 94, 87, 73, 8, 60, 70, 67…
#> $ dummy <lgl> TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE,…Coercion
Before utilising methods from qalytools we must first coerce our input data to the <EQ5D5L> object class via as_eq5d5l()
dat <- as_eq5d5l(
EQ5D5L_surveys,
mobility = "mobility",
self_care = "self_care",
usual = "usual",
pain = "pain",
anxiety = "anxiety",
respondentID = "respondentID",
surveyID = "surveyID",
vas = "vas",
)
dat
#> # EQ-5D-5L: 10,000 x 12
#> surve…¹ respo…² sex age mobil…³ self_…⁴ usual pain anxiety time_…⁵ vas
#> * <chr> <int> <chr> <int> <int> <int> <int> <int> <int> <dbl> <dbl>
#> 1 survey… 1 Fema… 25 2 1 2 1 1 30 96
#> 2 survey… 1 Fema… 25 1 4 5 5 3 60 22
#> 3 survey… 1 Fema… 25 1 5 2 1 3 90 61
#> 4 survey… 1 Fema… 25 1 2 5 2 2 120 65
#> 5 survey… 1 Fema… 25 1 2 3 3 1 150 78
#> 6 survey… 1 Fema… 25 2 1 1 2 1 180 88
#> 7 survey… 1 Fema… 25 2 1 2 1 1 210 85
#> 8 survey… 1 Fema… 25 1 1 2 2 1 240 95
#> 9 survey… 1 Fema… 25 1 3 2 1 1 270 94
#> 10 survey… 1 Fema… 25 2 1 2 1 1 300 87
#> # … with 9,990 more rows, 1 more variable: dummy <lgl>, and abbreviated
#> # variable names ¹surveyID, ²respondentID, ³mobility, ⁴self_care, ⁵time_index
#> # ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable namesDescriptive methods
To obtain a quick overview of the data we can call summary(). By default this returns the output as a list of data frames, showing frequency counts and proportions, split by surveyID
head(summary(dat), n = 2)
#> $survey01
#> # A data frame: 6 × 6
#> value mobility self_care usual pain anxiety
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 666 608 401 320 574
#> 2 2 196 131 303 359 235
#> 3 3 114 205 166 287 150
#> 4 4 19 18 61 24 29
#> 5 5 5 38 69 10 12
#> 6 NA 0 0 0 0 0
#>
#> $survey02
#> # A data frame: 6 × 6
#> value mobility self_care usual pain anxiety
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 396 106 142 26 47
#> 2 2 248 169 177 71 92
#> 3 3 194 214 160 93 103
#> 4 4 83 277 237 364 381
#> 5 5 79 234 284 446 377
#> 6 NA 0 0 0 0 0Alternatively we can set the parameter tidy to TRUE and obtain the summary data in a “tidy” (long) format.
summary(dat, tidy = TRUE)
#> # A data frame: 300 × 4
#> surveyID dimension value count
#> <chr> <fct> <dbl> <dbl>
#> 1 survey01 mobility 1 666
#> 2 survey01 mobility 2 196
#> 3 survey01 mobility 3 114
#> 4 survey01 mobility 4 19
#> 5 survey01 mobility 5 5
#> 6 survey01 mobility NA 0
#> 7 survey01 self_care 1 608
#> 8 survey01 self_care 2 131
#> 9 survey01 self_care 3 205
#> 10 survey01 self_care 4 18
#> # … with 290 more rows
#> # ℹ Use `print(n = ...)` to see more rowsIt can be useful to look at the percentage of responses, across dimensions, that are not equal to one. This can be done with the calculate_limitation() function
limitation <- calculate_limitation(dat)
limitation
#> # A data frame: 50 × 3
#> surveyID dimension without_limitation
#> <chr> <fct> <dbl>
#> 1 survey01 mobility 0.666
#> 2 survey02 mobility 0.396
#> 3 survey03 mobility 0.647
#> 4 survey04 mobility 0.637
#> 5 survey05 mobility 0.613
#> 6 survey06 mobility 0.637
#> 7 survey07 mobility 0.645
#> 8 survey08 mobility 0.651
#> 9 survey09 mobility 0.642
#> 10 survey10 mobility 0.671
#> # … with 40 more rows
#> # ℹ Use `print(n = ...)` to see more rows
ggplot(limitation, aes(x = surveyID, y = without_limitation, group = dimension)) +
geom_line(aes(colour = dimension)) +
theme_light() +
scale_y_continuous(n.breaks = 10, expand = c(0.005, 0.005), limits = c(0, 1)) +
scale_x_discrete() +
scale_fill_discrete(name = "dimension") +
ylab("Without limitation")
Calculating utility values
To calculate the utility values we call calculate_utility() on our EQ5D object with additional arguments specifying the countries and type we are interested in. calculate_utility() will return the desired values split by country, type, respondentID and surveyID. If you prefer to augment the utility values on top of the input object we also provide the add_utility() function.
We can obtain a list of compatible value sets (across countries and type) by passing our EQ5D object directly to the available_valuesets() function or by passing a comparable string.
For EQ5D3L inputs the type can be:
“TTO”, the time trade-off valuation technique;
“VAS”, the visual analogue scale valuation technique;
“RCW”, a reverse crosswalk conversion to EQ5D5L values; or
“DSU”, the NICE Decision Support Unit’s model that allows mappings on to EQ5D5L values accounting for both age and sex.
For EQ5D5L inputs this can be:
“VT”, value sets generated via a EuroQol standardised valuation study protocol;
“CW”, a crosswalk conversion EQ5D3L values; or
“DSU”, the NICE Decision Support Unit’s model that allows mappings on to EQ5D5L values accounting for both age and sex.
Note that available_valuesets()is a convenience wrapper around the eq5d::valuesets() function. It will return a data frame with columns representing the EQ5D version, the value set country and the value set type.
vs <- available_valuesets(dat)
head(vs, 10)
#> # A data frame: 10 × 3
#> Version Type Country
#> * <chr> <chr> <chr>
#> 1 EQ-5D-5L VT Belgium
#> 2 EQ-5D-5L VT Canada
#> 3 EQ-5D-5L VT China
#> 4 EQ-5D-5L VT Denmark
#> 5 EQ-5D-5L VT Egypt
#> 6 EQ-5D-5L VT England
#> 7 EQ-5D-5L VT Ethiopia
#> 8 EQ-5D-5L VT France
#> 9 EQ-5D-5L VT Germany
#> 10 EQ-5D-5L VT HongKong
# comparable character values will also give the same result
identical(vs, available_valuesets("eq5d5l"))
#> [1] TRUE
# cross walk comparison
vs <- filter(vs, Country %in% c("England", "UK"), Type %in% c("VT", "CW"))
util <- calculate_utility(dat, type = vs$Type, country = vs$Country)
# plot the results
util$fill <- paste(util$.utility_country, util$.utility_type)
ggplot(util, aes(x = surveyID, y = .value, fill = fill)) +
geom_boxplot(lwd = 1, outlier.shape = 4) +
stat_summary(
mapping = aes(group = .utility_country),
fun = mean,
geom = "point",
position = position_dodge(width = 0.75),
shape = 21,
color = "black",
fill = "white"
) +
theme_light() +
geom_hline(yintercept = 0, linetype = "longdash", size = 0.6, color = "grey30") +
scale_y_continuous(n.breaks = 10, expand = c(0.005, 0.005)) +
scale_x_discrete(name = "value set")
For further details about the available value sets, consult the documentation of the wrapped eq5d package via vignette(topic = "eq5d", package = "eq5d").
QALY calculations
Quality of life years can be calculated directly from utility values. By default, two different metrics are provided. Firstly, a “raw” value which is simply the scaled area under the utility curve and, secondly, a value which represents the loss from full health.
qalys <-
dat |>
add_utility(type = "VT", country = c("Denmark", "France")) |>
calculate_qalys(time_index = "time_index")
filter(qalys, .qaly=="raw")
#> # A data frame: 2,000 × 5
#> respondentID .utility_type .utility_country .qaly .value
#> <int> <chr> <chr> <chr> <dbl>
#> 1 1 VT Denmark raw 0.544
#> 2 1 VT France raw 0.573
#> 3 2 VT Denmark raw 0.487
#> 4 2 VT France raw 0.503
#> 5 3 VT Denmark raw 0.530
#> 6 3 VT France raw 0.577
#> 7 4 VT Denmark raw 0.276
#> 8 4 VT France raw 0.430
#> 9 5 VT Denmark raw 0.506
#> 10 5 VT France raw 0.575
#> # … with 1,990 more rows
#> # ℹ Use `print(n = ...)` to see more rows
filter(qalys, .qaly=="loss_vs_fullhealth")
#> # A data frame: 2,000 × 5
#> respondentID .utility_type .utility_country .qaly .value
#> <int> <chr> <chr> <chr> <dbl>
#> 1 1 VT Denmark loss_vs_fullhealth 0.195
#> 2 1 VT France loss_vs_fullhealth 0.166
#> 3 2 VT Denmark loss_vs_fullhealth 0.252
#> 4 2 VT France loss_vs_fullhealth 0.237
#> 5 3 VT Denmark loss_vs_fullhealth 0.210
#> 6 3 VT France loss_vs_fullhealth 0.162
#> 7 4 VT Denmark loss_vs_fullhealth 0.463
#> 8 4 VT France loss_vs_fullhealth 0.310
#> 9 5 VT Denmark loss_vs_fullhealth 0.233
#> 10 5 VT France loss_vs_fullhealth 0.164
#> # … with 1,990 more rows
#> # ℹ Use `print(n = ...)` to see more rowscalculate_qalys() also allows us to calculate the loss from a specified baseline in one of two ways. Firstly, a character string baseline_survey argument can be placed which matches a survey present in the utility data. If the argument is passed as such then the utility values from the specified survey are used to calculate the loss. Note that the survey is still included in the raw, unadjusted calculation, prior to the calculation of loss.
# Reload the example data and combine with some baseline measurements
data("EQ5D5L_surveys")
dat <-
EQ5D5L_surveys |>
as_eq5d5l(
mobility = "mobility",
self_care = "self_care",
usual = "usual",
pain = "pain",
anxiety = "anxiety",
respondentID = "respondentID",
surveyID = "surveyID",
vas = "vas"
)
add_utility(dat, type = "VT", country = c("Denmark", "France")) |>
calculate_qalys(baseline_survey = "survey01", time_index = "time_index") |>
filter(.qaly == "loss_vs_baseline")
#> # A data frame: 2,000 × 5
#> respondentID .utility_country .utility_type .qaly .value
#> <int> <chr> <chr> <chr> <dbl>
#> 1 1 Denmark VT loss_vs_baseline 0.141
#> 2 1 France VT loss_vs_baseline 0.114
#> 3 2 Denmark VT loss_vs_baseline 0.0892
#> 4 2 France VT loss_vs_baseline -0.00370
#> 5 3 Denmark VT loss_vs_baseline 0.0965
#> 6 3 France VT loss_vs_baseline 0.106
#> 7 4 Denmark VT loss_vs_baseline 0.175
#> 8 4 France VT loss_vs_baseline 0.0450
#> 9 5 Denmark VT loss_vs_baseline 0.142
#> 10 5 France VT loss_vs_baseline 0.0915
#> # … with 1,990 more rows
#> # ℹ Use `print(n = ...)` to see more rowsAlternatively the baseline_survey argument can be specified as a data frame with a column corresponding to the respondentID and another representing the associated utility. Optionally columns corresponding to the utility country and utility type can be included to allow more granular comparisons. Note that for this specification of baseline, it is not included in the unadjusted, raw, calculation.
split_dat <- split(dat, dat$surveyID=="survey01")
surveys <- split_dat[[1]]
baseline <- split_dat[[2]]
utility_dat <- add_utility(
surveys,
type = "VT",
country = c("Denmark", "France")
)
baseline_utility <-
baseline |>
add_utility(type = "VT", country = c("Denmark", "France")) |>
select(respondentID,.utility_country,.utility_type,.value)
calculate_qalys(utility_dat, baseline_survey = baseline_utility, time_index = "time_index") |>
filter(.qaly == "loss_vs_baseline")
#> # A data frame: 2,000 × 5
#> respondentID .utility_country .utility_type .qaly .value
#> <int> <chr> <chr> <chr> <dbl>
#> 1 1 Denmark VT loss_vs_baseline 0.101
#> 2 1 France VT loss_vs_baseline 0.0796
#> 3 2 Denmark VT loss_vs_baseline 0.0556
#> 4 2 France VT loss_vs_baseline -0.0341
#> 5 3 Denmark VT loss_vs_baseline 0.0672
#> 6 3 France VT loss_vs_baseline 0.0788
#> 7 4 Denmark VT loss_vs_baseline 0.145
#> 8 4 France VT loss_vs_baseline 0.0182
#> 9 5 Denmark VT loss_vs_baseline 0.0959
#> 10 5 France VT loss_vs_baseline 0.0553
#> # … with 1,990 more rows
#> # ℹ Use `print(n = ...)` to see more rowsOther functionality
Paretian Classification of Health Change (PCHC)
data("eq5d3l_example")
dat <- as_eq5d3l(
eq5d3l_example,
respondentID = "respondentID",
surveyID = "surveyID",
mobility = "MO",
self_care = "SC",
usual = "UA",
pain = "PD",
anxiety = "AD",
vas = "vas",
drop = FALSE # ensure we do not drop additional columns
)
grp1 <- filter(dat, Group == "Group1")
grp2 <- filter(dat, Group == "Group2")
calculate_pchc(grp1, grp2)
#> # A data frame: 5 × 3
#> Change Number Percent
#> <chr> <dbl> <dbl>
#> 1 No change 14 14
#> 2 Improve 59 59
#> 3 Worsen 14 14
#> 4 Mixed change 13 13
#> 5 No problems 0 0
calculate_pchc(grp1, grp2, by.dimension = TRUE)
#> # A data frame: 20 × 4
#> .Dimension Change Number Percent
#> <chr> <chr> <dbl> <dbl>
#> 1 MO No change 31 31
#> 2 MO Improve 26 26
#> 3 MO Worsen 7 7
#> 4 MO No problems 36 36
#> 5 SC No change 21 21
#> 6 SC Improve 31 31
#> 7 SC Worsen 7 7
#> 8 SC No problems 41 41
#> 9 UA No change 33 33
#> 10 UA Improve 43 43
#> 11 UA Worsen 6 6
#> 12 UA No problems 18 18
#> 13 PD No change 59 59
#> 14 PD Improve 34 34
#> 15 PD Worsen 3 3
#> 16 PD No problems 4 4
#> 17 AD No change 14 14
#> 18 AD Improve 22 22
#> 19 AD Worsen 12 12
#> 20 AD No problems 52 52