---
title: "Introduction"
date: "`r format(as.Date('2026-06-30'), '%d %B, %Y')`"
author:
  - name: "Samvel B. Gasparyan"
    affiliation: https://gasparyan.co/
output:
  rmarkdown::html_document:
    theme: "darkly"
    highlight: "zenburn"
    toc: true
    toc_float: true
link-citations: true
bibliography: REFERENCES.bib
biblio-style: science
vignette: >
  %\VignetteIndexEntry{Introduction}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style>
a[role="doc-biblioref"] {
  color: green !important;
}

a[role="doc-biblioref"]:visited {
  color: green !important;
}

a[role="doc-biblioref"]:hover {
  color: lightgreen !important;
  text-decoration: underline;
}
</style>

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
R <- function() knitr::include_graphics("Rlogo.png", dpi = 5000)
```

## `hce` package intro

### Background

```{r echo=FALSE, out.width = '33%'}
knitr::include_graphics("hex-hce.png")
```


The package supports the simulation and analysis of **univariate hierarchical composite endpoints (HCEs)**. The rationale for constructing a univariate endpoint, in contrast to the multivariate version of Generalized Pairwise Comparisons, which typically yields summary statistics rather than a single endpoint, is described in Gasparyan, Koch, and Brunner [@gasparyan2026univariate]. The framework is broad: univariate HCEs can be derived from different outcome types, including time-to-event, continuous, and binary outcomes. However, it applies only to fixed follow-up designs without early dropout.

In this setting, several approaches can be used to derive a patient-level ordinal outcome, which is the essence of a univariate HCE. Gasparyan, Koch, and Brunner [@gasparyan2026univariate] derived the distribution for the setting with multiple time-to-event outcomes and a single continuous (or ordinal) outcome assessed at the end of follow-up; in that construction, the event times also contribute to the HCE. Accordingly, our definition of an HCE focuses on this setting. Although the theoretical distribution is not required for estimation, understanding the distribution remains important for interpreting summary statistics in greater detail and for avoiding the **Condorcet non-transitivity paradox** [@de2014essai].

These endpoints have been implemented in clinical trials across multiple therapeutic areas. For example, see @gasparyan2022design for an implementation in a COVID-19 setting and for practical considerations in constructing hierarchical composite endpoints. For chronic kidney disease (**CKD**) applications, see [@khce1; @khce2; @khce3].

The primary analysis method is **win odds**. Other win statistics, including **win ratio** and **net benefit** [@dong2023win], are also implemented when there is no censoring. Win odds relies on the **DeLong-DeLong-Clarke-Pearson** formula [@delong1988comparing] for the variance of the win proportion and is based on the **Brunner-Munzel** test [@brunner2000nonparametric]. We also include the **Brunner-Konietschke** version [@brunner2025unbiased] of the **Bamber** estimator [@bamber1975area] for the variance of the win proportion. In addition, the package provides **Wilson-type, range-preserving confidence intervals** for the win proportion and win odds, following Schüürhuis, Konietschke, and Brunner [@schuurhuis2025new].

The power and sample size formulas cover several classes of alternatives: `"shifted"`, `"ordered"`, and `"max"` for the maximum value of the standard deviation. All formulas are derived from Bamber [@bamber1975area]. By default, **Noether's formula** [@noether1987sample] is used for shifted distributions. Additional discussion of power calculations for shifted distributions is available in [@gasparyan2021power; @gasparyan2022comments].

For reviews of HCE design in clinical trials, see [@gasparyan2022design], [@gasparyanhierarchical], and [@khce1]. The **Basic Data Structure (BDS)**, which follows **Analysis Data Model (ADaM)** [@CDISC] principles for hierarchical composite endpoints, is described in @gasp2024bds. For visualization, the **maraca plot** [@karpefors2023maraca] can be used.

Stratified and adjusted win odds are calculated using **randomization-based covariate adjustment theory** [@koch1998issues]; for a review, see [@gasparyan2021adjusted].

Thresholds for pairwise comparisons involving the continuous component of the HCE are implemented according to the theory in Gasparyan, Koch, and Brunner [@gasparyan2026univariate].

All implementations are **rank-based**. For a methodological review, see [@brunner2018rank].


### Setup

Load the package `hce` and check the version

```{r eval = TRUE}
library(hce)
packageVersion("hce")
```

For citing the package, run `citation("hce")` [@hce].

### Contents

List the functions and the datasets in the package

```{r}
ls("package:hce")
```

The package contains the following **Datasets**: use `data(package = "hce")` for the list of all datasets included in the package.

```{r eval=FALSE}
HCE1; HCE2; HCE3; HCE4
COVID19; COVID19b; COVID19plus
ADET; ADLB; ADSL; KHCE
```

- **Simulated datasets**: `HCE1 - HCE4` contain two treatment groups and analysis values `AVAL` for a hierarchical composite endpoint.

- **COVID-19 datasets**: `COVID19`, `COVID19b`, and `COVID19plus` contain COVID-19 ordinal scale outcomes [@beigel2020remdesivir; @kalil2021baricitinib].

- **Synthetic kidney disease datasets**: `ADET` contains event-time data, `ADLB` contains laboratory data, and `ADSL` contains subject-level baseline characteristics. Based on these data, `KHCE` provides the derived kidney hierarchical composite endpoint for the same patients [@khce2].

## In Brief

### Calculations

```{r eval=FALSE}
calcWO()
calcWINS()
summaryWO()
regWO()
stratWO()
IWP()
```

The main functions for calculating the so-called win statistics—**win probability**, **win odds**, **win ratio**, and **net benefit**—are `calcWO()` for win odds and `calcWINS()` for all win statistics, including win odds. To summarize wins, losses, and ties, use `summaryWO()`, which also reports the win probability, win odds, and win ratio, together with the standard errors of the win probability and win odds. By default, all of these functions use the **DeLong-DeLong-Clarke-Pearson** formula [@delong1988comparing] to estimate the standard error of the win probability and win odds.

In `calcWINS()`, the `SE_WP_Type` argument can be set to `unbiased` to use the **Brunner-Konietschke** version [@brunner2025unbiased] of the **Bamber** estimator [@bamber1975area] for the variance of the win proportion, or to `biased` to use the **DeLong-DeLong-Clarke-Pearson** variance estimator, which is the same method used in `calcWO()` and may be biased in small samples. In the same function, the `CI_WP_Type` argument can be set to `Wilson` to obtain **Wilson-type, range-preserving confidence intervals** for the win proportion and win odds, following Schüürhuis, Konietschke, and Brunner [@schuurhuis2025new]. The function `calcWINS()` also implements **Goodman-Kruskal's gamma** [@goodman1954measures; @goodman1963measures]. All reported statistics include confidence intervals.

The functions `regWO()` and `stratWO()` compute adjusted and stratified win odds, respectively, based on the **randomization-based covariate adjustment theory** developed in Koch et al. [@koch1998issues]; for a review, see Gasparyan et al. [@gasparyan2021adjusted].

In principle, all of these functions operate on ordinal outcomes when each patient contributes a single analysis value. Therefore, the input data should be at the patient level. Any patient-level dataset with a numeric analysis value can be used to calculate these statistics. In this setting, the use of win statistics is equivalent to ordinal analysis. Specifically, win odds corresponds to the classical **Mann-Whitney odds**, win ratio is the number of concordances divided by the number of discordances and can be derived from **Goodman-Kruskal's gamma**, and net benefit corresponds to **Somers' D C/R**. Therefore, any software that computes these statistics should produce the same results. For example, although the implementation here is fully rank-based, the estimate of **Goodman-Kruskal's gamma** and its confidence interval agree with `DescTools::GoodmanKruskalGamma()` [@Desc].

All calculations are rank-based and therefore computationally efficient. The basic building blocks are the **individual win proportions**, defined for each patient as the proportion of wins plus one-half of the proportion of ties. The function `IWP()` computes these individual win proportions for each patient in the dataset and returns an updated dataset with a new column. Any win-odds-related calculation can be based on these quantities [@gasparyan2021adjusted].

### Univariate hierarchical composite endpoints

```{r eval=FALSE}
hce()
as_hce()
```


The idea behind **univariate hierarchical composite endpoints** [@gasparyan2026univariate] is to derive a single ordinal outcome from multiple outcomes of different types, such as time-to-event, continuous, and binary outcomes, by prioritizing the most clinically important outcome for each patient. The resulting ordinal outcome can then be used to calculate win statistics.

To support this, the package provides helper functions for deriving `hce` objects from different types of outcomes. In particular, `hce()` and `as_hce()` check the structure of the input data and attempt to convert it to an `hce` object. Once this conversion is complete, that is, once the univariate HCE has been derived, the subsequent theory and analysis are the same as for any other ordinal outcome.

### Thresholds

```{r eval=FALSE}
deltaWO()
```


The function `deltaWO()` calculates threshold-adjusted win odds [@gasparyan2026univariate]. 

### Simulations

```{r eval=FALSE}
simORD()
simHCE()
simTTE()
simKHCE()
```


### Power and sample size


```{r eval=FALSE}
powerWO()
sizeWO()
sizeWR()
minWO()
propWINS()
```

The functions `powerWO()`, `sizeWO()`, and `minWO()` provide tools for calculating **power**, **sample size**, and the **minimum detectable treatment effect** for win odds under different classes of alternatives: `"shifted"`, `"ordered"`, and `"max"`, where `"max"` refers to the maximum value of the standard deviation. All formulas are based on Bamber [@bamber1975area]. For shifted distributions, the default is **Noether's formula** [@noether1987sample]. See also [@gasparyan2021power; @gasparyan2022comments] for further discussion of shifted distributions. The function `sizeWR()` provides sample size calculations for the **win ratio** [@yu2022sample].

The function `propWINS()` is a convenient tool for deriving the proportions of **wins**, **losses**, and **ties** for each treatment group from the win odds and win ratio.

### Ordinal dominance graph and the maraca plot

```{r eval=FALSE}
plot.hce()
```

A plot method for `hce` objects (created by `as_hce()`) to provide the ordinal dominance graph [@bamber1975area].

### Extras

```{r eval=FALSE}
rweibullGF()
dGLL()
pGLL()
qGLL()
rGLL()
hGLL()
HGLL()
```

The **generalized log-logistic (GLL)** distributions arise as the marginal distributions of the **Weibull-Gamma frailty model**. For comparison, one can use `rweibullGF()`, which implements the random frailty formulation, and compare its output with the generalized log-logistic distribution functions: density `dGLL()`, distribution function `pGLL()`, quantile function `qGLL()`, random generation `rGLL()`, hazard function `hGLL()`, and cumulative hazard function `HGLL()`.

These distributions form the basis for constructing two-outcome hierarchical composite endpoints with death and hospitalization in the generalized illness-death model implemented in `simTTE()`.

### Helpers

```{r eval=FALSE}
print.hce_results()
plot.hce_results()
as_formulae.formula()
```


To facilitate visualization and review of sample size and power calculations, `print()` and `plot()` methods are implemented for `hce_results` objects generated by the functions `powerWO()`, `sizeWO()`, and `minWO()`.

The function `as_formulae()` is an internal helper that decomposes formula objects into treatment, response, covariate, and grouping variables using the standardized names `TRTP`, `AVAL`, `COVAR`, and `GROUP`, respectively.

## `hce` Objects


### `hce()` Function

`hce` objects can be constructed using the helper function `hce()`, which has the following arguments:

```{r}
args("hce")
```
 
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':

```{r}
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)
```

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()`:

```{r}
head(dat)
```

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):

```{r}
unique(dat[, c("GROUP", "PARAMN")])
```

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:

```{r}
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")])
```

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:


```{r}
data(HCE1)
unique(HCE1$GROUP)
```

Therefore, we can construct the following object using the `hce()` function:

```{r}
HCE <- hce(GROUP = factor(HCE1$GROUP, levels = c("TTE1", "TTE2", "TTE3", "TTE4", "C")), 
           TRTP = HCE1$TRTP, AVAL0 = HCE1$AVAL0, PADY = 1080)
class(HCE)
head(HCE)
```

The object's class is `adhce`, which inherits from `hce`. This design indicates that the object includes additional columns.

### Create an `hce` Object from a Data Frame

Consider the dataset `HCE1`, which is part of the package `hce`:

```{r}
data(HCE1, package = "hce")
class(HCE1)
head(HCE1)
```

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.

```{r}
dat1 <- as_hce(HCE2)
str(dat1)
```


### Simulate `hce` Objects Using `simHCE()`

To simulate values from a hierarchical composite endpoint, we use the function `simHCE()`, which has the following arguments:

```{r}
args("simHCE")
```

- 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.

```{r}
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)
```

```{r}
class(dat3)
head(dat3)
```




## Generics 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:
```{r}
methods(class = "hce")
```

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.

```{r}
HCE <- hce(GROUP = factor(HCE3$GROUP, levels = c("TTE1", "TTE2", "TTE3", "TTE4", "C")), 
           TRTP = HCE3$TRTP, AVAL0 = HCE3$AVAL0, PADY = 1080)
calcWO(HCE)
calcWINS(HCE)  
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:

```{r}
methods(class = "adhce")
```

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:

```{r}
summaryWO(HCE)
```

An important addition to `summaryWO.adhce()` is the `cumsummary_by_GROUP` output, which reports cumulative wins, losses, and ties as outcomes are added sequentially according to the priority order:

```{r}
res0 <- summaryWO(HCE, ref = "P")
res <- res0$cumsummary_by_GROUP
barplot(PROP ~ WINS + GROUPN, data = res, 
col = c("darkgreen", "darkred", "darkblue"), 
xlab = "Proportions", xlim = c(0, 1), 
ylab = "Cumulative components by prioritization", 
legend.text = unique(res$WINS), beside = TRUE, horiz = TRUE)
grid()
```


## References