hce
package introThe 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).
Load the package hce
and check the version
For citing the package, run citation("hce")
(Samvel B. Gasparyan 2025).
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:
data(package = "hce")
for the list of all
datasets included in the package.Simulated datasets - HCE1 - HCE4
that contain two
treatment groups and analysis values AVAL
of a hierarchical
composite endpoint.
The datasets COVID19
, COVID19b
,
COVID19plus
of 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 dataset
KHCE
for the same patients (Heerspink et al. 2023).
hce
objects: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 an
hce
object for the kidney HCE (Heerspink et al. 2023).
The function simORD()
simulates ordinal outcomes by
categorization of beta distributions.
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).gamma
) and their confidence intervals calculation
function:calcWINS()
(see Samvel B.
Gasparyan et al. (2023)). The argument SE_WP_Type
by
default is set to biased
since it produces
DeLong-DeLong-Clarke-Pearson variance (the same as in
calcWO()
) which is biased for small samples. The
SE_WP_Type
can be changed to unbiased
to
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()
.
regWO()
(see Samvel B.
Gasparyan, Folkvaljon, et al. (2021); Samvel B. Gasparyan et al. (2023)).stratWO()
(see Samvel B.
Gasparyan, Folkvaljon, et al. (2021); Samvel B. Gasparyan et al. (2023)).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_results
objects,
generated by the functions
powerWO(), sizeWO(), minWO()
.
A plot method for hce
objects (created by
as_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
Objectshce()
FunctionThe 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 the following arguments:
We 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 hce
structure assumes that only one event per
patient is present for the analysis, meaning that the resulting
hce
object created by the hce()
function is a
patient-level dataset. The function hce()
does not select
the clinically most important event of the patient but requires it to be
already done when calling it.
The argument TRTP
should 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 PARAMN AVAL
#> 1 Control I 1 1
#> 2 Active I 1 1
#> 3 Control I 1 1
#> 4 Active I 1 1
#> 5 Active I 1 1
#> 6 Control I 1 1
We see that the dataset has a very specific structure. The column
PARAMN
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 GROUP
argument can be provided as a factor with the levels specifying the
necessary order:
set.seed(2022)
n <- 100
GROUP = rep(x = c("I", "II", "III"), each = 100)
GROUP <- factor(GROUP, levels = c("III", "II", "I"))
dat <- hce(GROUP = GROUP,
TRTP = sample(x = c("Active", "Control"), size = n*3, replace = TRUE))
unique(dat[, c("GROUP", "PARAMN")])
#> GROUP PARAMN
#> 1 I 3
#> 101 II 2
#> 201 III 1
This 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 object using the
hce()
function:
HCE <- hce(GROUP = factor(HCE1$GROUP, levels = c("TTE1", "TTE2", "TTE3", "TTE4", "C")),
TRTP = HCE1$TRTP, AVAL0 = HCE1$AVAL0, PADY = 1080)
class(HCE)
#> [1] "adhce" "hce" "data.frame"
head(HCE)
#> TRTP GROUP AVAL0 PARAMN AVAL PADY GROUPN
#> 1 A C -2.21 5 5446.32 1080 5400
#> 2 P C 17.06 5 5465.59 1080 5400
#> 3 A TTE4 966.00 4 5286.00 1080 4320
#> 4 P TTE4 352.00 4 4672.00 1080 4320
#> 5 A C -12.45 5 5436.08 1080 5400
#> 6 P C 44.62 5 5493.15 1080 5400
The object’s class is adhce
, which inherits from
hce
. This design indicates that the object includes
additional columns.
hce
Object from a Data FrameConsider 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 1080
This 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 ...
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)
#> NULL
The vector arguments TTE_A
and TTE_P
specify 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 default
shape = 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 argument
pat
. The function simulates event times in days, and the
yeardays = 360
argument 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 fixedfy
argument 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 6
hce
ObjectsAs 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 code
The 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 = factor(HCE3$GROUP, levels = c("TTE1", "TTE2", "TTE3", "TTE4", "C")),
TRTP = HCE3$TRTP, AVAL0 = HCE3$AVAL0, PADY = 1080)
calcWO(HCE)
#> WO LCL UCL SE WOnull alpha Pvalue WP
#> 1 1.332829 1.15078 1.543678 0.07493202 1 0.05 0.0001014229 0.571336
#> LCL_WP UCL_WP SE_WP SD_WP N
#> 1 0.5353673 0.6073047 0.01835169 0.5803314 1000
calcWINS(HCE)
#> $summary
#> WIN LOSS TIE TOTAL Pties
#> 1 142824 107156 20 250000 8e-05
#>
#> $WP
#> WP LCL UCL Pvalue
#> 1 0.571336 0.5353673 0.6073047 0.0001014229
#>
#> $NetBenefit
#> NetBenefit LCL UCL Pvalue
#> 1 0.142672 0.0707347 0.2146093 0.0001014229
#>
#> $WO
#> WO LCL UCL Pvalue
#> 1 1.332829 1.15078 1.543678 0.0001014229
#>
#> $WR1
#> WR LCL1 UCL1 Pvalue1
#> 1 1.332861 1.150793 1.543733 0.000125974
#>
#> $WR2
#> WR LCL2 UCL2 Pvalue2
#> 1 1.332861 1.155092 1.537987 8.351682e-05
#>
#> $gamma
#> gamma LCL UCL Pvalue
#> 1 0.1426834 0.07074057 0.2146263 0.000101418
#>
#> $SE
#> WP_SE NetBenefit_SE logWR_SE1 logWR_SE2 gamma_SE
#> 1 0.01835169 0.03670338 0.07493804 0.07303552 0.03670621
#>
#> $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)
To check the generic functions available for the adhce
class by running the following command:
The main difference between summaryWO.hce()
and
summaryWO.adhce()
is that summaryWO.adhce()
also provides a group-wise summary, in addition to the overall
summary:
summaryWO(HCE)
#> $summary
#> TRTP WIN LOSS TIE TOTAL WR WO
#> 1 A 142824 107156 20 250000 1.3328605 1.3328294
#> 2 P 107156 142824 20 250000 0.7502661 0.7502835
#>
#> $summary_by_GROUP
#> TRTP GROUP WIN LOSS TIE TOTAL
#> 1 A TTE1 2473 39525 2 42000
#> 2 P TTE1 2985 29513 2 32500
#> 3 A TTE2 6847 22648 5 29500
#> 4 P TTE2 13857 45138 5 59000
#> 5 A TTE3 13251 16745 4 30000
#> 6 P TTE3 13356 25140 4 38500
#> 7 A TTE4 9056 6443 1 15500
#> 8 P TTE4 15383 19616 1 35000
#> 9 A C 111197 21795 8 133000
#> 10 P C 61575 23417 8 85000
#>
#> $WO
#> WO SE WP SE_WP
#> 1 1.332829 0.07493202 0.571336 0.01835169