hce package intro
Background
The purpose of the package is to simulate and analyze hierarchical composite endpoints (HCEs). The primary analysis method is win odds, but other win statistics (win ratio, net benefit)(Dong et al. 2023) are also implemented, provided there is no censoring. Win odds uses the DeLong-DeLong-Clarke-Pearson fromula (DeLong, DeLong, and Clarke-Pearson 1988) for the variance of the win proportion and is based on the Brunner-Munzel test (Brunner and Munzel 2000).
The power and sample size formulas consider different alternative classes (“shifted”, “ordered”, and “max” for the maximum value of the standard deviation). All formulas are derived from Bamber (1975). By default, Noether’s formula (Noether 1987) is used for shifted distributions. For more information on the power calculations for shifted distributions see also Samvel B. Gasparyan, Kowalewski, et al. (2021), Samvel B. Gasparyan, Kowalewski, and Koch (2022).
For a review of designing HCEs in clinical trials, see Samvel B. Gasparyan et al. (2022), also Samvel B. Gasparyan et al. (2023) and Little et al. (2023). The Basic Data Structure (BDS), which conforms to the Analysis Data Model (ADaM) (CDISC 2009) principles for hierarchical composite endpoints, is detailed in Samvel B. Gasparyan et al. (2024). To visualize HCEs, the maraca plot (Karpefors, Lindholm, and Gasparyan 2023) can be utilized.
The stratified and adjusted win odds are calculated based on the randomization-based covariate adjustment theory developed in Koch et al. (1998) (for a review, see Samvel B. Gasparyan, Folkvaljon, et al. (2021)).
All implementations are rank-based. For the review of methodology see (Brunner, Bathke, and Konietschke 2018).
Setup
Load the package hce and check the version
For citing the package, run citation("hce") (Samvel B. Gasparyan 2025).
Contents
List the functions and the datasets in the package
ls("package:hce")
#> [1] "ADET" "ADLB" "ADSL" "COVID19" "COVID19b"
#> [6] "COVID19plus" "HCE1" "HCE2" "HCE3" "HCE4"
#> [11] "IWP" "KHCE" "as_hce" "calcWINS" "calcWO"
#> [16] "hce" "minWO" "powerWO" "propWINS" "regWO"
#> [21] "simHCE" "simKHCE" "simORD" "sizeWO" "sizeWR"
#> [26] "stratWO" "summaryWO"In brief, the package contains the following:
- Datasets: use
data(package = "hce")for the list of all datasets included in the package.
Simulated datasets -
HCE1 - HCE4that contain two treatment groups and analysis valuesAVALof a hierarchical composite endpoint.The datasets
COVID19,COVID19b,COVID19plusof COVID-19 ordinal scale outcomes (Beigel et al. 2020; Kalil et al. 2021).The datasets
ADET(event-time),ADLB(laboratory),ADSL(subject-level baseline characteristics) for kidney events and their timing, kidney-related laboratory measurements of eGFR (estimated glomerular filtration rate), and, based on this, the derived kidney hierarchical composite endpoint datasetKHCEfor the same patients (Heerspink et al. 2023).
- Functions to create
hceobjects:
hce(), as_hce(), simHCE()(see Samvel B. Gasparyan et al. (2024)).The function
simKHCE()simulates continuous eGFR (estimated Glomerular Filtration Rate) values over time and derives anhceobject for the kidney HCE (Heerspink et al. 2023).The function
simORD()simulates ordinal outcomes by categorization of beta distributions.
- Win odds (win ratio with ties) calculation generic functions for hierarchical composite endpoints:
calcWO(), summaryWO()(see Samvel B. Gasparyan, Kowalewski, et al. (2021)). Win odds uses the DeLong-DeLong-Clarke-Pearson formula (DeLong, DeLong, and Clarke-Pearson 1988) for the variance of the win proportion and is based on the Brunner-Munzel test (Brunner and Munzel 2000).
- Win statistics (win odds, win ratio, net benefit, Goodman Kruskal’s
gamma) and their confidence intervals calculation function:
calcWINS()(see Samvel B. Gasparyan et al. (2023)). The argumentSE_WP_Typeby default is set tobiasedsince it produces DeLong-DeLong-Clarke-Pearson variance (the same as incalcWO()) which is biased for small samples. TheSE_WP_Typecan be changed tounbiasedto produce Brunner-Konietschke version (Brunner and Konietschke 2025) of the Bamber (Bamber 1975) estimator for the variance of win proportion.Back-calculation of number of wins, losses, and ties given the win odds and win ratio using the function
propWINS().
- Adjusted win odds calculation for a single, numeric covariate:
regWO()(see Samvel B. Gasparyan, Folkvaljon, et al. (2021); Samvel B. Gasparyan et al. (2023)).
- Stratified win odds with a possible adjustment after stratification with a single, numeric covariate:
stratWO()(see Samvel B. Gasparyan, Folkvaljon, et al. (2021); Samvel B. Gasparyan et al. (2023)).
- Power, sample size, and minimum detectable win odds calculation functions:
powerWO(), sizeWO(), minWO()provide tools for the win odds power, sample size, and minimum detectable treatment effect calculation for different alternative classes (“shifted”, “ordered”, and “max” for the maximum value of the standard deviation) based on the win odds. All formulas are from Bamber (1975). By default, uses Noether’s formula (Noether 1987) for shifted distributions. For shifted distributions see also Samvel B. Gasparyan, Kowalewski, et al. (2021), Samvel B. Gasparyan, Kowalewski, and Koch (2022).Win ratio sample size calculation formula
sizeWR()(Yu and Ganju 2022).
Print and plot methods for
hce_resultsobjects, generated by the functionspowerWO(), sizeWO(), minWO().A plot method for
hceobjects (created byas_hce()) to provide the ordinal dominance graph (Bamber 1975).The function
IWP()to calculate patient-level individual win proportions (Samvel B. Gasparyan, Folkvaljon, et al. 2021; Samvel B. Gasparyan, Kowalewski, et al. 2021).
hce Objects
hce() Function
The main objects in the package are called hce objects,
which are data frames with a specific structure matching the design of
hierarchical composite endpoints (HCE). These are complex
endpoints that combine different clinical events into a composite, using
a hierarchy to prioritize the clinically most important event for a
patient. These endpoints are implemented in clinical trials across
different therapeutic areas. See, for example, Samvel B. Gasparyan et al. (2022) for an
implementation in a COVID-19 setting and some practical considerations
for constructing hierarchical composite endpoints. For the Chronic
Kidney Disease (CKD) outcomes see (Little et al.
2023; Heerspink et al. 2023; Kondo et al. 2024).
HCEs are ordinal endpoints that can be thought of as having ‘greater’, ‘less’, or ‘equal’ defined for them, but without a definition of how much greater or less. In this sense, the ordinal outcomes can be represented as numeric vectors as long as numeric operations (e.g., sum or division) are not performed on them.
hce objects can be constructed using the helper function
hce(), which has three arguments:
args("hce")
#> function (GROUP = character(), TRTP = character(), AVAL0 = 0,
#> ORD = sort(unique(GROUP)))
#> NULLWe see that the required arguments are GROUP, which
specifies the clinically most important event of a patient to be
included in the analysis, and TRTP, which specifies the
(planned) treatment group of a patient (exactly two treatment groups
should be present). Note that:
The
hcestructure assumes that only one event per patient is present for the analysis, meaning that the resultinghceobject created by thehce()function is a patient-level dataset. The functionhce()does not select the clinically most important event of the patient but requires it to be already done when calling it.The argument
TRTPshould have exactly two levels.
Consider the following example of ordinal outcomes ‘I’, ‘II’, and ‘III’:
set.seed(2022)
n <- 100
dat <- hce(GROUP = rep(x = c("I", "II", "III"), each = 100),
TRTP = sample(x = c("Active", "Control"), size = n*3, replace = TRUE))
class(dat)
#> [1] "hce" "data.frame"This dataset has the appropriate structure of hce
objects, but its class inherits from an object of class
data.frame. This means that all functions available for
data frames can be applied to hce objects, for example, the
function head():
head(dat)
#> TRTP GROUP GROUPN AVAL AVAL0 ORD
#> 1 Control I 0 0 0 1
#> 2 Active I 0 0 0 1
#> 3 Control I 0 0 0 1
#> 4 Active I 0 0 0 1
#> 5 Active I 0 0 0 1
#> 6 Control I 0 0 0 1We see that the dataset has a very specific structure. The column
GROUPN shows how the function hce() generated
the order of given events (it uses usual alphabetic order for the unique
values in the GROUP column to determine the clinical
importance of events):
In the class hce, higher values for the ordering mean
clinically less important events. For example, death, which is the most
important event, should always get the lowest ordinal value. If there is
a need to specify the order of outcomes, then the argument
ORD can be used:
set.seed(2022)
n <- 100
dat <- hce(GROUP = rep(x = c("I", "II", "III"), each = 100),
TRTP = sample(x = c("A", "P"), size = n*3, replace = TRUE), ORD = c("III", "II", "I"))
unique(dat[, c("GROUP", "GROUPN")])
#> GROUP GROUPN
#> 1 I 2
#> 101 II 1
#> 201 III 0This means that the clinically most important event is ‘III’ instead
of ‘I’. The argument AVAL0 is meant to help in cases where
we want to introduce sub-ordering within each GROUP
category. For example, if two events in the group ‘I’ can be compared
based on other parameters, then the AVAL0 argument can be
specified to take that into account.
Below we use the built-in data frame HCE1 to construct
an hce object. Before specifying the order of events, it is
a good idea to check what are the unique events included in the
GROUP column:
Therefore, we can construct the following hce
object:
HCE <- hce(GROUP = HCE1$GROUP, TRTP = HCE1$TRTP, AVAL0 = HCE1$AVAL0,
ORD = c("TTE1", "TTE2", "TTE3", "TTE4", "C"))
class(HCE)
#> [1] "hce" "data.frame"
head(HCE)
#> TRTP GROUP GROUPN AVAL AVAL0 ORD
#> 1 A C 40000 39997.79 -2.21 10000
#> 2 P C 40000 40017.06 17.06 10000
#> 3 A TTE4 30000 30966.00 966.00 10000
#> 4 P TTE4 30000 30352.00 352.00 10000
#> 5 A C 40000 39987.55 -12.45 10000
#> 6 P C 40000 40044.62 44.62 10000Create an hce Object from a Data Frame
Consider the dataset HCE1, which is part of the package
hce:
data(HCE1, package = "hce")
class(HCE1)
#> [1] "data.frame"
head(HCE1)
#> ID TRTP GROUP GROUPN AVALT AVAL0 AVAL PADY
#> 1 1 A C 5400 1080 -2.21 5446.32 1080
#> 2 2 P C 5400 1080 17.06 5465.59 1080
#> 3 3 A TTE4 4320 966 966.00 5286.00 1080
#> 4 4 P TTE4 4320 352 352.00 4672.00 1080
#> 5 5 A C 5400 1080 -12.45 5436.08 1080
#> 6 6 P C 5400 1080 44.62 5493.15 1080This dataset has the appropriate structure of hce
objects, but its class is data.frame. A generic way of
coercing data structures to an hce object is to use the
function as_hce(). This function performs checks (using an
internal validator function) and creates an hce object from
the given data structure (using an internal constructor function). If
coercion is not possible, it will throw an error explaining the
issue.
dat1 <- as_hce(HCE2)
str(dat1)
#> Classes 'adhce', 'hce' and 'data.frame': 1000 obs. of 9 variables:
#> $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ TRTP : chr "A" "P" "A" "P" ...
#> $ GROUP : Factor w/ 5 levels "TTE1","TTE2",..: 1 5 5 1 2 3 5 2 1 5 ...
#> $ GROUPN: num 1080 5400 5400 1080 2160 3240 5400 2160 1080 5400 ...
#> $ AVALT : num 120 1080 1080 577 782 985 1080 293 313 1080 ...
#> $ AVAL0 : num 120 3.35 22.8 577 782 ...
#> $ AVAL : num 1200 5461 5480 1657 2942 ...
#> $ PADY : num 1080 1080 1080 1080 1080 1080 1080 1080 1080 1080 ...
#> $ PARAMN: num 1 5 5 1 2 3 5 2 1 5 ...Simulate hce Objects Using simHCE()
To simulate values from a hierarchical composite endpoint, we use the
function simHCE(), which has the following arguments:
args("simHCE")
#> function (n, n0 = n, TTE_A, TTE_P, CM_A, CM_P, CSD_A = 1, CSD_P = CSD_A,
#> fixedfy = 1, yeardays = 360, pat = 100, shape = 1, theta = 1,
#> logC = FALSE, seed = NULL, dec = 2, all_data = FALSE)
#> NULLThe vector arguments
TTE_AandTTE_Pspecify the event rates per year for time-to-event outcomes in the active and control groups, respectively. The function assumes a Weibull survival function with the same shape parameter for simulating all time-to-event outcomes in both treatment groups (by defaultshape = 1, which assumes an exponential survival function). These two vectors should have the same length, which indicates the number of time-to-event outcomes.By default, the event rates are presented per 100 patient-years (
pat = 100), which can be changed using the argumentpat. The function simulates event times in days, and theyeardays = 360argument can be used to change the number of days in a year (e.g., 365 or 365.25).The function simulates events during a fixed follow-up period only, and the
fixedfyargument can be used to change the length of the follow-up (in years).The function simulates the continuous outcome from a normal (default) or log-normal (if
logC = TRUE) distribution with given means and standard deviations for two treatment groups.
Rates_A <- c(1.72, 1.74, 0.58, 1.5, 1)
Rates_P <- c(2.47, 2.24, 2.9, 4, 6)
dat3 <- simHCE(n = 2500, n0 = 1500, TTE_A = Rates_A,
TTE_P = Rates_P,
CM_A = -3, CM_P = -6,
CSD_A = 16, CSD_P = 15,
fixedfy = 3, seed = 2023)class(dat3)
#> [1] "adhce" "hce" "data.frame"
head(dat3)
#> ID TRTP GROUP GROUPN AVALT AVAL0 AVAL seed PADY PARAMN
#> 1 1 A C 6480 1080 -9.32 6526.85 2023 1080 6
#> 2 2 A C 6480 1080 0.11 6536.28 2023 1080 6
#> 3 3 A C 6480 1080 27.06 6563.23 2023 1080 6
#> 4 4 A C 6480 1080 -4.78 6531.39 2023 1080 6
#> 5 5 A TTE2 2160 390 390.00 2550.00 2023 1080 2
#> 6 6 A C 6480 1080 4.56 6540.73 2023 1080 6Generics for hce Objects
As we see, the function simHCE() creates an object of
type hce, which inherits from the built-in class
data.frame. We can check all implemented methods for this
new class as follows:
methods(class = "hce")
#> [1] calcWINS calcWO plot summaryWO
#> see '?methods' for accessing help and source codeThe function calcWO() calculates the win odds and its
confidence interval, while summaryWO() provides a more
detailed calculation of win odds, including the number of wins, losses,
and ties by GROUP categories.
HCE <- hce(GROUP = HCE3$GROUP, TRTP = HCE3$TRTP,
ORD = c("TTE1", "TTE2", "TTE3", "TTE4", "C"))
calcWO(HCE)
#> WO LCL UCL SE WOnull alpha Pvalue WP
#> 1 1.333635 1.160574 1.532502 0.07091636 1 0.05 3.852522e-05 0.571484
#> LCL_WP UCL_WP SE_WP SD_WP N
#> 1 0.5374459 0.6055221 0.01736671 0.5491836 1000
summaryWO(HCE)
#> $summary
#> TRTP WIN LOSS TIE TOTAL WR WO
#> 1 A 110655 74913 64432 250000 1.4771135 1.3336352
#> 2 P 74913 110655 64432 250000 0.6769961 0.7498303
#>
#> $summary_by_GROUP
#> TRTP GROUP WIN LOSS TIE TOTAL
#> 1 A C 87780 0 45220 133000
#> 2 P C 39780 0 45220 85000
#> 3 A TTE1 0 36540 5460 42000
#> 4 P TTE1 0 27040 5460 32500
#> 5 A TTE2 3835 18703 6962 29500
#> 6 P TTE2 9912 42126 6962 59000
#> 7 A TTE3 10980 14400 4620 30000
#> 8 P TTE3 11011 22869 4620 38500
#> 9 A TTE4 8060 5270 2170 15500
#> 10 P TTE4 14210 18620 2170 35000
#>
#> $WO
#> WO SE WP SE_WP
#> 1 1.333635 0.07091636 0.571484 0.01736671
calcWINS(HCE)
#> $summary
#> WIN LOSS TIE TOTAL Pties
#> 1 110655 74913 64432 250000 0.257728
#>
#> $WP
#> WP LCL UCL Pvalue
#> 1 0.571484 0.5374459 0.6055221 3.852522e-05
#>
#> $NetBenefit
#> NetBenefit LCL UCL Pvalue
#> 1 0.142968 0.07489175 0.2110443 3.852522e-05
#>
#> $WO
#> WO LCL UCL Pvalue
#> 1 1.333635 1.160574 1.532502 3.852522e-05
#>
#> $WR1
#> WR LCL1 UCL1 Pvalue1
#> 1 1.477113 1.222982 1.784053 5.129812e-05
#>
#> $WR2
#> WR LCL2 UCL2 Pvalue2
#> 1 1.477113 1.226016 1.779638 4.069634e-05
#>
#> $gamma
#> gamma LCL UCL Pvalue
#> 1 0.1926086 0.1017117 0.2835056 3.279693e-05
#>
#> $SE
#> WP_SE NetBenefit_SE logWR_SE1 logWR_SE2 gamma_SE
#> 1 0.01736671 0.03473342 0.09632721 0.09506293 0.04637683
#>
#> $ref
#> [1] "A vs P"
#>
#> $Input
#> alpha WOnull
#> 1 0.05 1
HCE$TRTP <- factor(HCE$TRTP, levels = c("P", "A"))
plot(HCE, fill = TRUE, col = "#865A4F", type = 'l', lwd = 2)
abline(a = 0, b = 1, lwd = 2, col = "#999999", lty = 2)