Covariate-adjusted win statistics in randomized clinical trials with ordinal outcomes

TL;DR

This paper introduces covariate-adjusted estimators for ordinal outcomes in clinical trials using win statistics, improving efficiency over unadjusted methods and providing practical tools for analysis.

Contribution

It develops propensity score weighting and augmented estimators for win parameters in ordinal outcomes, with theoretical properties and implementation in R.

Findings

Weighted estimators outperform unadjusted methods in simulations

Augmented estimators further improve efficiency in most cases

Methods are demonstrated with the ORCHID trial data

Abstract

Ordinal outcomes are common in clinical settings where they often represent increasing levels of disease progression or different levels of functional impairment. Such outcomes can characterize differences in meaningful patient health states that are directly relevant to clinical researchers and frequently represent composite outcomes that include absorbing states such as death. To compare different intervention strategies in clinical trials, the direct use of ordinal logistic regression models may not be ideal for analyzing ranked outcomes due to non-collapsibility, lack of estimation and clarity, or failure of the common underlying proportional odds assumption. In this article, we focus on representing the average treatment effect for ordinal outcomes via intrinsic pairwise outcome comparisons captured through win estimates, such as the win ratio and win difference. We first develop propensity score weighting estimators, including both inverse probability weighting (IPW) and overlap weighting (OW), tailored to estimating win parameters. Furthermore, we develop augmented weighting estimators that leverage an additional ordinal outcome regression to potentially improve efficiency over weighting alone. Leveraging the theory of U-statistics, we establish the asymptotic theory for all estimators, and derive closed-form variance estimators to support statistical inference. Through extensive simulations we demonstrate the enhanced efficiency of the weighted estimators over the unadjusted estimator, with the augmented weighting estimators showing a further improvement in efficiency except for extreme cases. Finally, we illustrate our proposed methods with the ORCHID trial, and implement our covariate adjustment methods in an R package winPSW to facilitate the practical implementation.

Equations

$${\tau }_1=\mathbb{P}\left\{Y_i(1)>Y_j(0)\right\},$$

$$\theta =\frac{{\tau }_1}{{\tau }_{-1}}=\frac{\mathbb{P}\left\{Y_i(1)>Y_j(0)\right\}}{\mathbb{P}\left\{Y_i(1)<Y_j(0)\right\}},\eta ={\tau }_1-{\tau }_{-1}=\mathbb{P}\left\{Y_i(1)>Y_j(0)\right\}-\mathbb{P}\left\{Y_i(1)<Y_j(0)\right\},$$

$${\hat{\tau }}_1^{\text{UNADJ}}=\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{Z}_i(1-Z_j)\mathbb{I}(Y_i>Y_j)}{{\sum }_{i=1}^n{\sum }_{j=1}^n{Z}_i(1-Z_j)}.$$

$${\hat{\tau }}_1^{\text{UNADJ}}=\frac{1}{n_1{n}_0}\sum _{i=1}^{n_1}\sum _{j=1}^{n_0}\frac{Z_i(1-Z_j)\mathbb{I}(Y_i>Y_j)}{{n_1}^{-1}{\sum }_{i=1}^{n_1}{Z}_i\times {n_0}^{-1}{\sum }_{j=1}^{n_0}(1-Z_j)}=\frac{1}{n_1{n}_0}\sum _{i=1}^{n_1}\sum _{j=1}^{n_0}{\omega }_{ij}\mathbb{I}(Y_i>Y_j),$$

$$\sqrt{n}\left({\hat{\tau }}_1^{\text{UNADJ}}-{\tau }_1\right)\stackrel{d}{\sim }N\left(0,{\sigma }_1^2\right),$$

$${\xi }_{10}$$

$${\xi }_{01}$$

$${\xi }_{10}^{12}$$

$${\xi }_{01}^{12}$$

$$\text{var}(\hat{\theta })=\text{var}({\hat{\tau }}_1^{\text{UNADJ}}/{\hat{\tau }}_{-1}^{\text{UNADJ}})\approx \frac{({\hat{\tau }}_1^{\text{UNADJ}}{)}^2}{({\hat{\tau }}_{-1}^{\text{UNADJ}}{)}^2}[\frac{\text{var}({\hat{\tau }}_1^{\text{UNADJ}})}{({\hat{\tau }}_1^{\text{UNADJ}}{)}^2}-2\frac{\text{cov}({\hat{\tau }}_1^{\text{UNADJ}},{\hat{\tau }}_{-1}^{\text{UNADJ}})}{{\hat{\tau }}_1^{\text{UNADJ}}\cdot {\hat{\tau }}_{-1}^{\text{UNADJ}}}+\frac{\text{var}({\hat{\tau }}_{-1}^{\text{UNADJ}})}{({\hat{\tau }}_{-1}^{\text{UNADJ}}{)}^2}].$$

$$\text{var}(\hat{\eta })=\text{var}({\hat{\tau }}_1^{\text{UNADJ}})+\text{var}({\hat{\tau }}_{-1}^{\text{UNADJ}})-2\text{cov}({\hat{\tau }}_1^{\text{UNADJ}},{\hat{\tau }}_{-1}^{\text{UNADJ}}).$$

$${\tau }_1^h=\frac{\mathbb{E}\left\{h(X_i,X_j)\mathbb{I}\left[Y_i(1)>Y_j(0)\right]\right\}}{\mathbb{E}\left[h(X_i,X_j)\right]}.$$

$${\hat{\tau }}_1^h=\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n\omega (X_i,X_j)Z_i(1-Z_j)\mathbb{I}(Y_i>Y_j)}{{\sum }_{i=1}^n{\sum }_{j=1}^n{Z}_i(1-Z_j)\omega (X_i,X_j)},$$

$${\hat{\tau }}_1^h\stackrel{p}{⟶}{\tau }_1^h,\text{as}\text{penalty}10000n\to \infty .$$

$${\hat{\tau }}_1^{\text{IPW}}=\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{Z}_i(1-Z_j)\mathbb{I}(Y_i>Y_j)/\left\{\hat{e}(X_i)[1-\hat{e}(X_j)]\right\}}{{\sum }_{i=1}^n{\sum }_{j=1}^n{Z}_i(1-Z_j)/\left\{\hat{e}(X_i)[1-\hat{e}(X_j)]\right\}}.$$

$${\hat{\tau }}_1^{\text{OW}}=\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{Z}_i(1-Z_j)[1-\hat{e}(X_i)]\hat{e}(X_j)\mathbb{I}(Y_i>Y_j)}{{\sum }_{i=1}^n{\sum }_{j=1}^n{Z}_i(1-Z_j)[1-\hat{e}(X_i)]\hat{e}(X_j)}.$$

$${\tau }^g=\int \int g(t,s){\nu }_1(dt){\nu }_0(ds),$$

$${\hat{\tau }}_1^h={\left(\genfrac{}{}{0ex}{}{n}{2}\right)}^{-1}\sum _{i=1}^n\sum _{j>i}^n{q}^h({\mathcal{O}}_i,{\mathcal{O}}_j;\hat{\beta }),$$

$$q^h({\mathcal{O}}_i,{\mathcal{O}}_j;\hat{\beta })=\frac{n(n-1)\left\{Z_i(1-Z_j)\mathbb{I}(Y_i>Y_j)\omega (X_i,X_j)+Z_j(1-Z_i)\mathbb{I}(Y_i<Y_j)\omega (X_j,X_i)\right\}}{2{\sum }_{i=1}^n{\sum }_{j=1}^n{Z}_i(1-Z_j)\omega (X_i,X_j)}.$$

$${\varphi }^h(\mathcal{O})=2\{q_1^h(\mathcal{O};\beta )-{\tau }_1^h\}+\frac{\partial }{\partial {\beta }^{\top }}(\mathbb{E}\times \mathbb{E})q^h(\mathcal{O},\mathcal{O};\beta ){\tilde{l}}_{\beta }(\mathcal{O};\beta ),$$

$${\chi }^h(\mathcal{O})$$

$$n^{-2}\sum _{i=1}^n{\left[2\left\{{\hat{q}}_1^h({\mathcal{O}}_i;{\hat{\beta }}_n)-{\hat{\tau }}_1^h\right\}+{\tilde{l}}_{\beta }^{\top }(Z_i,X_i;\hat{\beta })A_1^h\right]}^{\otimes 2},$$

$$\widehat{\text{var}}\left({\hat{\tau }}_1^{\text{UNADJ}}\right)=n^{-2}\sum _{i=1}^n{\left[2\left\{{\hat{q}}_1^{\text{UNADJ}}({\mathcal{O}}_i)-{\hat{\tau }}_1^{\text{UNADJ}}\right\}\right]}^{\otimes 2},$$

$${\hat{\tau }}_{-1}^h=\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n\omega (X_i,X_j)Z_i(1-Z_j)\mathbb{I}(Y_i<Y_j)}{{\sum }_{i=1}^n{\sum }_{j=1}^n\omega (X_i,X_j)Z_i(1-Z_j)},$$

$$\mu (X_i,X_j;\alpha )=\mathbb{E}[g\{Y_i(1),Y_j(0)\}|X_i,X_j;\alpha ]=\iint g(t,s){\nu }_1(dt|X_i;\alpha ){\nu }_0(ds|X_j;\alpha ).$$

$${\hat{\tau }}_1^{\text{MAO}}={\left(\genfrac{}{}{0ex}{}{n}{2}\right)}^{-1}\sum _{i=1}^n\sum _{j>i}^n{r}^{\text{MAO}}({\mathcal{O}}_i,{\mathcal{O}}_j;\hat{\beta },\hat{\alpha }),$$

$$r^{\text{MAO}}({\mathcal{O}}_i,{\mathcal{O}}_j;\hat{\beta },\hat{\alpha })$$

$${\hat{\tau }}_1^{\text{MAO}}$$

$$\mu (X_i,X_j;\alpha )$$

$$\text{logit}(P(Y\le l|X))={\alpha }_{l0}-{\xi }^{\top }X\iff P(Y\le l|X)=\frac{\exp ({\alpha }_{l0}-{\xi }^{\top }X)}{1+\exp ({\alpha }_{l0}-{\xi }^{\top }X)}.$$

$$p_l(X;\alpha )=P(Y=l)=\{\begin{array}{ll}P(Y\le 1|X)=P(Y=1|X),& \text{if }Y=l=1\\ P(Y\le l|X)-P(Y\le l-1|X)& \text{if }1<Y=l<L\\ 1-P(Y\le L-1|X),& \text{if }Y=l=L.\end{array}$$

Full text

Covariate-adjusted win statistics in randomized clinical trials with ordinal outcomes

Zhiqiang Cao1, Scott Zuo2, Mary Ryan Baumann3,4, Kendra Plourde5,6, Patrick Heagerty7,8,

Guangyu Tong5,9,10,11 and Fan Li5,9,10,11,∗

1School of Artificial Intelligence, Shenzhen Technology University, Guangdong, P. R. China

2Department of Preventive Medicine, Feinberg School of Medicine, Northwestern University, Illinois, U.S.A

3Department of Biostatistics and Medical Informatics, University of Wisconsin - Madison, Wisconsin, U.S.A

4Department of Population Health Sciences, University of Wisconsin - Madison, Wisconsin, U.S.A

5Department of Biostatistics, Yale University School of Public Health, Connecticut, U.S.A

6Yale Center for Analytical Sciences, Yale University School of Public Health, Connecticut, U.S.A

7Department of Biostatistics, University of Washington, Washington, U.S.A

8Clinical Learning, Evidence, and Research (CLEAR) Center, University of Washington, Washington, U.S.A

9Section of Cardiovascular Medicine, Department of Internal Medicine, Yale School of Medicine, Connecticut, U.S.A

10Cardiovascular Medicine Analytics Center, Yale School of Medicine, Connecticut, U.S.A

11Center for Methods in Implementation and Prevention Science, Yale School of Public Health, Connecticut, U.S.A

$*$fan.f.li@yale.edu

Abstract

Ordinal outcomes are common in clinical settings where they often represent increasing levels of disease progression or different levels of functional impairment. In this article, we focus on representing the average treatment effect for ordinal outcomes via intrinsic pairwise outcome comparisons captured through win estimands, such as the win ratio and win difference. Recognizing the value of baseline covariate adjustment toward enhanced precision, we first develop propensity score weighting estimators, including both inverse probability weighting (IPW) and overlap weighting (OW), tailored to estimating win estimands. Furthermore, we develop augmented weighting estimators that leverage an additional ordinal outcome regression to potentially improve efficiency over weighting alone. Leveraging the theory of U-statistics, we establish the asymptotic theory for all estimators, and derive closed-form variance estimators to support statistical inference. We also prove that all of the covariate-adjusted estimators do not compromise consistency for the target estimand even when the associated working models are incorrectly specified; hence these covariate-adjusted estimators are model-robust. Through simulations we demonstrate the enhanced efficiency of the weighted estimators over the unadjusted estimator, with the augmented weighting estimators showing a further improvement in efficiency except for extreme cases. Finally, we illustrate our proposed methods with the ORCHID trial, and implement our covariate adjustment methods in an R package winPSW.

Key words: Win ratio, Net benefit, inverse probability weighting, overlap weighting, Augmented weighting, U-statistics

1 Introduction

Ordinal outcomes are common in randomized clinical trials (RCTs) and are often comprised of multiple, monotonically ordered categories that do not necessarily align with a quantifiable distance (Selman et al., 2023). For example, the FLU-IVIG trial (Peterson et al., 2017) considered an ordinal outcome constructed with six mutually exclusive categories: (1) death; (2) intensive care unit (ICU) hospitalization (in ICU); (3) non-ICU hospitalization, requiring supplemental oxygen; (4) non-ICU hospitalization, not requiring supplemental oxygen; (5) discharged from the hospital but unable to resume normal activities; and (6) discharged from the hospital with resumption of normal activities. The categories were defined to delineate clear improvement and worsening in patient status and to yield a sufficient spread of the data for showing benefit due to intravenous hyperimmune immunoglobulin. Ordinal outcomes are useful in clinical settings, as specific disease states can represent meaningfully distinct categories that may be of clinical importance to researchers. However, the statistical considerations for analyzing ordinal outcomes can often be more complicated, especially when it comes to defining an assumption-lean target parameter to quantify the treatment effect. Conventional methods are often based on parametric ordinal regression models (e.g., proportional odds model) with an odds ratio target parameter (Uddin et al., 2024); however, if the proportional odds assumption fails to hold, the estimated odds ratio may be subject to bias with ambiguous interpretation. Instead, alternative summary estimands/measures such as the win ratio (Pocock et al., 2012; Mao, 2018; Klok and Siegerink, 2023), win odds (WO) (Dong et al., 2020b) — a modification of WR to include ties, and win difference (WD, or net benefit) (Buyse, 2010) have gained traction in recent years. Evans et al. (2020) showed that the win ratio can more optimally inform patient treatment, enhance the understanding of the totality of intervention effects on patients, and potentially provide efficiencies over standard analyses in a randomized clinical trial.

In this article, we focus on win statistics-based analyses of ordinal outcomes in randomized trials, and study covariate adjustment strategies to improve statistical efficiency. In RCTs, pre-specified covariate adjustment is a critical strategy for improving the precision of the treatment effect estimator as well as a crucial step in trials using randomization strata. Regulatory agencies such as the Food and Drug Administration (FDA) (U.S. Food and Drug Administration, 2023a) and the European Medicines Agency (EMA) (ICH, 2021) have explicitly recommended incorporating baseline covariates to enhance statistical efficiency in treatment comparisons. A rich literature has examined covariate adjustment methods for estimating the average treatment effects in randomized trials. For example, Benkeser et al. (2021) demonstrated that adjusting for baseline covariates in the analysis of binary, ordinal, and time-to-event outcomes can markedly reduce the required sample size by improving statistical precision. Shen et al. (2014) and Williamson et al. (2014) showed that inverse probability weighting (IPW) estimators, when applied with prognostic covariates to estimate known treatment propensities under randomization, can reduce the asymptotic variance compared to the unadjusted treatment effect estimator. Zeng et al. (2021) proposed overlap weighting (OW) for covariate adjustment in randomized trials, showing that OW estimators are asymptotically equivalent to IPW and ANCOVA estimators under 1:1 randomization but can yield lower finite-sample variance than IPW. Wang et al. (2023) established the consistency and asymptotic normality of a general class of M-estimators that integrate covariate adjustment and stratified randomization, showing substantial gains in efficiency and power. In contrast, far less attention has been devoted to covariate adjustment for pairwise comparison estimands, such as the win ratio. Notable among the sparse contributions, Vermeulen et al. (2015) extended the Mann-Whitney U test using semiparametric efficiency theory to exploit covariate information, leading to improved power and enabling exact permutation-based inference. Zhang et al. (2019) investigated several estimators—including regression-based, IPW, and augmented inverse probability weighting (AIPW) approaches—for Mann-Whitney-type causal effects in both clinical trials and observational studies. More recently, Wang et al. (2025) proposed inverse probability of treatment weighting (IPTW)-adjusted win ratio estimators for composite time-to-event outcomes, demonstrating improved power under non-proportional hazards. Buyse et al. (2025) briefly discussed covariate adjustment methods including IPW and AIPW for generalized pairwise comparison estimands in observational studies. While these previous efforts illustrate the potential of covariate adjustment (e.g., IPW) in pairwise comparisons, there is a lack of formal theoretical justification for the asymptotic efficiency gain of IPW in this context. Moreover, a more systematic investigation of a broader class of propensity score weighting or augmented weighting estimators for covariate-adjusted win analysis, particularly for ordinal outcomes, remains unavailable in the existing literature.

To address this gap, we study a general class of propensity score weighting methods for baseline covariate adjustment when using win statistics to analyze RCTs with ordinal outcomes. The key contributions of this work are threefold. First, we propose a class of covariate-adjusted weighting estimators for win statistics based on pairwise comparisons, incorporating both IPW and OW as two special cases, thereby extending the framework of Zeng et al. (2021) to a broader class of estimands. Second, we discuss augmented weighting estimators that combine any propensity score balancing weights with an outcome model to improve efficiency. This generalizes the AIPW formulation of Mao et al. (2019) to a wider range of balancing weighting schemes. Third, we establish the asymptotic properties of these estimators under randomized assignment by leveraging the theory of M-estimation and U-statistics. We further derive a set of closed-form variance estimators to enable computationally efficient inference without the need for resampling. Through this effort, a key theoretical insight is that the influence function of a balancing weighting estimator can be expressed as the influence function of the unadjusted estimator minus its projection onto the nuisance tangent space, for estimating the pairwise comparison estimands. This decomposition allows us to formally characterize the asymptotic variance reduction achieved through covariate adjustment for the win difference and win ratio, and to justify the use of propensity score weighting when applying win statistics in randomized trials. Through extensive simulations, we demonstrate that the proposed covariate-adjusted estimators, particularly augmented weighting estimators under either IPW or OW, substantially improve the efficiency of the unadjusted estimators for win estimands in randomized trials. We also implement the proposed covariate-adjusted win statistics in an R package winPSW to facilitate practical implementation.

The rest of this article is organized as follows. In Section 2, we introduce notation and estimands. In Section 3, we develop covariate-adjusted IPW and OW win estimators and derive their asymptotic properties. We study augmented weighting estimators that leverage an additional outcome regression to further improve efficiency in Section 4. Simulation studies are carried out to examine the finite-sample performance of the proposed estimators relative to the unadjusted estimator in Section 5. Finally, we illustrate our methods by reanalyzing the ordinal outcome data from the Outcomes Related to COVID-19 Treated With Hydroxychloroquine Among Inpatients With Symptomatic Disease (ORCHID) trial in Section 6, followed by a discussion in Section 7.

2 Win estimands and unadjusted estimators

2.1 Notation and estimands

We consider a two-arm RCT with a total of $n$ patients. Let $Z$ be the binary treatment assignment, where $Z=1$ for treatment and $Z=0$ for control. Let $n_{1}=\sum_{i=1}^{n}Z_{i}$ and $n_{0}=n-n_{1}$, denote the sample size in each arm. Suppose the outcome $Y$ is ordinal with sample space $\mathcal{Y}=\{1,2,\dots,L\}$, where $L$ is a positive finite integer. We assume that larger values of $Y$ are more favorable. Under the potential outcomes framework, each individual has a pair of potential outcomes $\{Y_{i}(1),Y_{i}(0)\}$, and we denote the observed outcome as $Y_{i}=Z_{i}Y_{i}(1)+(1-Z_{i})Y_{i}(0)$. Let $X=(X_{1},\dots,X_{p})^{\top}$ be a $p$-dimensional vector of baseline covariates collected for each patient. The observed data then consists of the tuple $\mathcal{O}_{i}=(Y_{i},Z_{i},X_{i}),\penalty 10000\ i=1,\dots,n$, which are $n$ independent and identically distributed realizations of $(Y,Z,X)$. With a comparison operator “$>$” or “$<$”, the win estimand is given by Wang and Pocock (2016):

[TABLE]

The above estimand has also been referred to by alternative terminology, including the Mann-Whitney effect, relative effect, and desirability of outcome ranking (DOOR) probability (Buyse et al., 2025); here we refer it as the win probability — that is, the probability that a randomly selected $Y_{i}(1)$ in the study population achieves a more favorable outcome than a randomly selected $Y_{j}(0)$. If there are no ties or ties are split evenly in both win and loss groups, a value greater than $1/2$ indicates that the treatment leads to more favorable outcomes stochastically. We can similarly define the loss estimand $\tau_{-1}=\mathbb{P}\left\{Y_{i}(1)<Y_{j}(0)\right\}$ and write the win ratio (WR) and win difference (WD, or net benefit) estimands, respectively, as:

[TABLE]

In what follows, we primarily discuss the estimation and inference of the win estimand and the associated covariate adjustment methods, before returning to the WR and WD estimands.

2.2 Unadjusted estimator

In a two-sample setting, Bebu and Lachin (2016) proposed an unbiased covariate-unadjusted estimator for the win estimand $\tau_{1}$ and loss estimand $\tau_{-1}$, and derived the large-sample results using multivariate multi-sample U-statistic. Below, we review and adapt their results under the potential outcomes framework. Under randomization, the treatment assignment $Z_{i}$ follows a Bernoulli random variable with probability $\pi=\mathbb{P}\left[Z_{i}=1|X_{i},Y_{i}(1),Y_{i}(0)\right]=\mathbb{P}(Z_{i}=1)\in(0,1)$. Under 1:1 randomization, $\pi=1/2$. Then the unadjusted win probability estimator is:

[TABLE]

Dividing both the numerator and denominator of (2.3) by $n^{2}$, we have $\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}Z_{i}(1-Z_{j})\mathbb{I}(Y_{i}>Y_{j})\stackrel{{\scriptstyle p}}{{\longrightarrow}}\mathbb{E}\left[Z_{i}(1-Z_{j})\mathbb{I}\left(Y_{i}(1)>Y_{j}(0)\right)\right]=\mathbb{E}\left[Z_{i}(1-Z_{j})\right]\mathbb{P}\left\{Y_{i}(1)>Y_{j}(0)\right\}$, and $\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}Z_{i}(1-Z_{j})\stackrel{{\scriptstyle p}}{{\longrightarrow}}\mathbb{E}\left[Z_{i}(1-Z_{j})\right]$, where $\mathbb{I}(\cdot)$ is an indicator function; then it is easy to show that $\widehat{\tau}_{1}^{\text{UNADJ}}$ is a consistent estimator of $\tau_{1}$.

We note that the index $i$ and $j$ in (2.3) are defined relative to the whole randomized sample. In practice, individual $i$ or $j$ can only belong to one of the treatment and control groups. If we treat $Y_{i}(1)$ and $Y_{j}(0)$ as coming from two groups with sample sizes equal to $n_{1}$ and $n_{0}$, respectively (similar to the setup in Bebu and Lachin (2016)), then $\widehat{\tau}_{1}^{\text{UNADJ}}$ is a multivariate multi-sample U-statistics:

[TABLE]

where the weight is given by $\omega_{ij}=\frac{Z_{i}(1-Z_{j})}{{n_{1}}^{-1}\sum_{i=1}^{n_{1}}Z_{i}\times{n_{0}}^{-1}\sum_{j=1}^{n_{0}}(1-Z_{j})}$ and (2.4) can be regarded as the weighted form of $U_{1}$ in (2.2) of Bebu and Lachin (2016). Then, according to Bebu and Lachin (2016), we have

[TABLE]

where $\sigma_{1}^{2}=\frac{n}{n_{1}}\xi_{10}+\frac{n}{n_{0}}\xi_{01}$, $n=n_{1}+n_{0}$ and

[TABLE]

Furthermore, by replacing “$>$” with “$<$” in (2.4), we obtain the corresponding unadjusted estimator for the loss probability estimand $\tau_{-1}$, denoted as $\widehat{\tau}_{-1}^{\text{UNADJ}}$, with similar large sample properties. Finally, the asymptotic covariance between $\widehat{\tau}_{1}^{\text{UNADJ}}$ and $\widehat{\tau}_{-1}^{\text{UNADJ}}$ is given by $\sigma_{12}=\frac{n}{n_{1}}\xi^{12}_{10}+\frac{n}{n_{0}}\xi^{12}_{01}$, where

[TABLE]

Once $\sigma_{12}$ is estimated, the covariance between $\widehat{\tau}_{1}^{\text{UNADJ}}$ and $\widehat{\tau}_{-1}^{\text{UNADJ}}$ can be computed as ${\text{cov}}\left(\widehat{\tau}_{1}^{\text{UNADJ}},\widehat{\tau}_{-1}^{\text{UNADJ}}\right)=\frac{1}{n}\widehat{\sigma}_{12}$. By applying a Taylor expansion, the variance of the estimated WR, $\widehat{\theta}$, based on unadjusted estimator (2.4) can be approximated as

[TABLE]

Note that the variance of the corresponding estimated WD, $\widehat{\eta}$, can be directly calculated as

[TABLE]

Hence, variance estimates of WR and WD from the unadjusted estimators can be obtained from (2.6) and (2.7), respectively. In Web Appendix A1-A2, we also present details about how to obtain consistent estimators for ${\text{var}}\left(\widehat{\tau}_{1}^{\text{UNADJ}}\right)$, ${\text{var}}\left(\widehat{\tau}_{-1}^{\text{UNADJ}}\right)$, and ${\text{cov}}\left(\widehat{\tau}_{1}^{\text{UNADJ}},\widehat{\tau}_{-1}^{\text{UNADJ}}\right)$.

Remark** 1**

Our $\widehat{\tau}_{1}^{\text{UNADJ}}$ in (2.4) is constructed following the general principles in Bebu and Lachin (2016) but differs from that existing estimator due to the weight $\omega_{ij}$. For both the variance expression in Bebu and Lachin (2016) and (2.5), we constrain the comparisons in the variance estimation to strictly different individuals (Web Appendix A1-A2), i.e., for both probabilities $\mathbb{P}\left[Y_{i}(1)>Y_{j}(0)\penalty 10000\ \&\penalty 10000\ Y_{i}(1)<Y_{k}(0)\right]$ and $\mathbb{P}\left[Y_{i}(1)>Y_{j}(0)\penalty 10000\ \&\penalty 10000\ Y_{k}(1)<Y_{j}(0)\right]$, we need $j\neq k$ and $i\neq k$. Ozenne et al. (2021) showed that this can lead to bias in small samples and potentially negative variance estimates for $\hat{\tau}_{1}^{\text{UNADJ}}$. They suggested a bias-reduced variance estimator based on first and second order H-decomposition of U-statistics that allows comparisons in the variance estimation to also come from the same individual (i.e., in the first probability, $j$ can be equal to $k$; and in the second probability, $i$ can be equal to $k$), which has less bias and can only be positive. Following their approach, we also provided a bias-reduced variance estimator for $\hat{\tau}_{1}^{\text{UNADJ}}$ in Web Appendix A3. By comparing the expressions given in Web Appendix A1 and A3, one can observe that the key difference lies in estimating the two probabilities $\mathbb{P}\left[Y_{i}(1)>Y_{j}(0)\penalty 10000\ \&\penalty 10000\ Y_{i}(1)<Y_{k}(0)\right]$ and $\mathbb{P}\left[Y_{i}(1)>Y_{j}(0)\penalty 10000\ \&\penalty 10000\ Y_{k}(1)<Y_{j}(0)\right]$ (by using (2.5) versus the method suggested by Ozenne et al. (2021)). We will also examine their performance by simulations in Section 5.

3 Covariate-adjusted win statistics via propensity score weighting

3.1 The balancing weights

In this Section, we introduce propensity score weighting methods for covariate-adjusted estimation of WRs and WDs. We first generalize the win estimand in (2.1) to a class of weighted win estimands to motivate alternative weighting estimators. Specifically, we assume that the marginal density of $X$ (from the combined treatment and control groups) exists and is denoted by $f(X)$. Under the balancing weight framework (Li et al., 2018; Li and Li, 2019), we can represent the target population density by $f(X)h(X)/C$, where $h(X)$ is a pre-specified tilting function, and $C$ is a normalization constant such that the density is well-defined and integrates to unity. The weighted win estimand for a given target population characterized by the tilting function is

[TABLE]

Here $h(X_{i},X_{j})$ is a bivariate tilting function defined as the ratio of the covariate distribution density in the target population to that in the study sample. In principle, $h(X_{i},X_{j})$ is a symmetric function of the two individuals’ covariates. In practice and throughout what follows, we will specify $h$ as a product of univariate tilting functions as $h(X_{i},X_{j})=h^{*}(X_{i})h^{*}(X_{j})$, where $h^{*}$ is the more familiar individual tilting function. This factorization satisfies the required symmetry and is natural given $h^{*}$ can be defined as a function of the propensity score, $e(x)=\mathbb{P}(Z=1|X)$.

In observational studies, propensity scores are usually unknown and must be estimated from observed data. Therefore, different $h(\cdot)$ corresponds to different populations and estimands. However, under randomization, for any baseline covariate value $X$, we have $e(X)=\mathbb{P}(Z=1)=\pi$.Thus, as long as $h(X_{i},X_{j})$ depends on covariates only through the propensity score, different $h(\cdot)$ correspond to the same target population and the weighted win estimand reduces to the unweighted win estimand: $\tau_{1}^{h}=\tau_{1}$. This unique feature provides the basis for considering alternative weighting estimators that perform implicit covariate adjustment without changing the target of inference (Zeng et al., 2021).

Zeng et al. (2021) leveraged the class of balancing weights to perform covariate adjustment, with the weights defined as $\omega_{1}(X)=h(X)/e(X)$ and $\omega_{0}=h(X)/\left[1-e(X)\right]$ for the treatment and control arms, respectively. We generalize their approach to handle the bivariate tilting function $h(X_{i},X_{j})$. Specifically, we focus on the following Hájek-type estimator:

[TABLE]

where $\omega(X_{i},X_{j})=\frac{\widehat{h}(X_{i},X_{j})}{\widehat{e}(X_{i})[1-\widehat{e}(X_{j})]}$. By the weak law of large numbers, we can show that $\widehat{\tau}_{1}^{h}$ is consistent to the target $\tau_{1}^{h}$ defined in (3.1):

[TABLE]

The results for $\tau_{-1}$ are analogous and are omitted for brevity.

3.2 Covariate adjustment via inverse probability weighting and overlap weighting

To leverage information from the baseline covariates, we first propose the following IPW estimator:

[TABLE]

Note that (3.4) is the special case of (3.2) by letting $\widehat{h}(X_{i},X_{j})=1$ and $\omega(X_{i},X_{j})=\frac{1}{\widehat{e}(X_{i})[1-\widehat{e}(X_{j})]}$. In practice, the propensity score $e(X)$ is usually estimated by a logistic regression model. Furthermore, it has been established that IPW can perform covariate adjustment and reduces the impact of chance imbalance, hence improving the precision over the unadjusted estimator, even when estimating a known propensity score rather than using the true propensity score $\pi$ (Shen et al., 2014; Williamson et al., 2014). Mao (2018) explored an IPW method for a general class of causal estimands defined by a contrast function, which extends the familiar notion of average treatment effect. However, their IPW estimator is unstabilized Horvitz-Thompson estimator, and hence differs from the stabilized IPW estimator (3.4). Wang et al. (2025) developed three IPW estimators: IPW-average treatment effect (IPW-ATE), stabilized IPW-ATE (SIPW-ATE) and IPW-average treatment effect in the treated (IPW-ATT) for the win ratio estimand. Nevertheless, their stabilized IPW-ATE is not a Hájek-type estimator as (3.4) and they generally ignored the uncertainty for estimating the propensity score in variance estimation. More recently, Buyse et al. (2025) briefly discussed covariate adjustment strategies including an IPW estimator for generalized pairwise comparison estimands, but they did not demonstrate how to conduct inference, nor did they provide a variance estimator. In contrast, we focus on a Hájek-type IPW estimator (3.4), and will derive a consistent variance estimator accounting for the uncertainty in estimating the propensity score.

Conventional IPW assigns a weight of $1/e(X)$ for those in the treatment arm and $1/\left[1-e(X)\right]$ for those in the control arm. Beyond IPW, OW assigns a weight of $1-e(X)$ for treated patients and a weight of $e(X)$ for control patients, and obviates the need for inverse weighting. In observational studies, it has been shown that OW retains the advantages of IPW while offering better finite-sample properties (Mao et al., 2019; Li and Li, 2019; Li et al., 2019; Cheng et al., 2022; Cao et al., 2025). By choosing $\widehat{h}(X_{i},X_{j})=\widehat{e}(X_{i})\widehat{e}(X_{j})(1-\widehat{e}(X_{i}))(1-\widehat{e}(X_{j}))$ so that $\omega(X_{i},X_{j})=\widehat{e}(X_{j})[1-\widehat{e}(X_{i})]$ in (3.2), we propose the Hájek-type estimator with OW as:

[TABLE]

In observational studies, OW has been shown to have the exact balance property and the minimum variance property (Li et al., 2018). That is, when $e(X)$ is estimated from a logistic regression, the propensity weighted covariate means are identical between the treatment and control groups (exact balance); under regularity conditions, OW leads to the most efficient estimator among the family of balancing weighting schemes for moment estimators (minimum asymptotic variance). The latter property also holds for modeling survival outcomes (Cheng et al., 2022; Cao et al., 2025). For conventional average treatment effect estimands, Zeng et al. (2021) further proved that OW and IPW lead to the same asymptotic variance in RCTs. However, due to the exact balance property of OW (Li et al., 2018), Zeng et al. (2021) demonstrated that OW can effectively remove all chance imbalance when the data generating process follows a linear model with additive terms. Furthermore, their simulation studies have shown that OW is always at least as efficient as IPW in finite samples for covariate adjustment, and carries a particular efficiency advantage when the sample size is small. Due to these reasons, we also investigate the use of OW in the context of WR and WD estimands.

3.3 Asymptotic analysis of the covariate-adjusted win statistics

Building off the theory of Mao (2018), we now apply U-statistic theory to derive the asymptotic variance of IPW and OW estimators. This derivation serves two purposes: (1) to show that propensity score weighting generally reduces the asymptotic variance when estimating win estimand, and (2) to obtain consistent variance estimators for covariate-adjusted win estimators. For convenience, we denote the propensity score model $e(X)$ by $\mathcal{M}_{e}$; here, we consider $\mathcal{M}_{e}$ to be a logistic regression — that is, $e(X)=\mathbb{P}(Z=1|X;\beta)=(1+\exp(-X^{\top}\beta))^{-1}$, where $\beta=(\beta_{1},\dots,\beta_{p})^{\top}$ is a finite $p$-dimensional parameter. Let $\widehat{\beta}$ be the maximum likelihood estimator for $\beta$ under $\mathcal{M}_{e}$, then $\widehat{e}(X)=(1+\exp(-X^{\top}\widehat{\beta}))^{-1}$. Let $g:\mathcal{X}\times\mathcal{X}\rightarrow\mathcal{R}$ be an arbitrary contrast function for two potential outcomes, where $\mathcal{X}$ is the space in which the potential outcomes reside, and let $\nu_{z}$ denote the marginal distribution of $Y(z)$, for $z=0,1$. Then, the win estimand takes the general form of

[TABLE]

based on a kernel $g$. When $g(t,s)=\mathbb{I}(t>s)$, $\tau^{g}=\tau_{1}$; similarly, when $g(t,s)=\mathbb{I}(t<s)$, $\tau^{g}=\tau_{-1}$. For generality, we focus on the asymptotic analysis for a general $g$ in what follows.

Since both the IPW estimator $\widehat{\tau}_{1}^{\text{IPW}}$ defined in (3.4) and the OW estimator $\widehat{\tau}_{1}^{\text{OW}}$ defined in (3.5) are special cases of the weighting estimator $\widehat{\tau}_{1}^{h}$ defined in (3.2), we can rewrite the weighting estimator as

[TABLE]

where

[TABLE]

Then, for any $i\neq j$, we have $\mathbb{E}\left[q^{h}(\mathcal{O}_{i},\mathcal{O}_{j};\widehat{\beta})\right]=\tau_{1}^{h}$ under randomization. By the Uniform Law of Large Numbers for the U-process (Arcones and Giné, 1993; Van Der Vaart and Wellner, 1996) and consistency of $\widehat{\beta}$, we can establish the consistency of $\widehat{\tau}_{1}^{h}$.

In Theorem 1, we show that: (1) the asymptotic variance of weighting estimator $\widehat{\tau}_{1}^{h}$ can be estimated by the second moment estimator for an estimated influence function; (2) both IPW and OW estimators have smaller asymptotic variance than the unadjusted estimator $\tau_{1}^{\text{UNADJ}}$. The proof of Theorem 1 is provided in Web Appendix B.

Theory** 1**

*(Asymptotic results of propensity score weighting estimators)

(a) The asymptotic variance of $n^{1/2}(\widehat{\tau}_{1}^{h}-\tau_{1}^{h})$ is ${\text{var}}\left(\phi^{h}(\mathcal{O})\right)$, where*

[TABLE]

*is the influence function of $\widehat{\tau}_{1}^{h}$, $q_{1}^{h}(\mathcal{O};\beta)=\mathbb{E}q^{h}\left(\cdot,\mathcal{O};\beta\right)=\frac{1}{2}\left[\frac{Z}{e(X)}\int\mathbb{I}(Y>s)\nu_{0}(ds)+\frac{1-Z}{1-e(X)}\int\mathbb{I}(Y<t)\nu_{1}(dt)\right]$, and $\widetilde{l}_{\beta}\left(\mathcal{O};\beta\right)$ is the influence function for $\beta$ based on logistic regression.

(b) The influence function $\phi^{h}(\mathcal{O})$ in (3.7) is equal to $\chi^{h}(\mathcal{O})$ minus its projection on the tangent space of $\mathcal{M}_{e}$, where*

[TABLE]

*The detailed expressions of the tangent space of $\mathcal{M}_{e}$ and the projection are provided in Web Appendix B.

(c) The unadjusted estimator $\widehat{\tau}_{1}^{\text{UNADJ}}$ also belongs to the same class of balancing weighting estimators $\widehat{\tau}_{1}^{h}$ defined in (3.2), and the influence function of $\widehat{\tau}_{1}^{\text{UNADJ}}$ is equal to $\chi^{h}(\mathcal{O})$ with $e(X)=\pi$. So the asymptotic variance of unadjusted estimator $\tau_{1}^{\text{UNADJ}}$ is equal to ${\text{var}}\left(\chi^{h}(\mathcal{O})\right)$. ***

Theorem 1(a) is obtained from applying Lemma S1 of Mao (2018) directly. By applying the same technique used in Remark S1 of Mao (2018) and the semi-parametric efficiency theory (Tsiatis, 2006), we can derive the result in Theorem 1(b). Note that if we let $\widehat{h}(X_{i},X_{j})=\widehat{e}(X_{i})[1-\widehat{e}(X_{j})]$ and $\omega(X_{i},X_{j})=1$, then $\widehat{\tau}_{1}^{h}$ defined in (3.2) is equivalent to $\tau_{1}^{\text{UNADJ}}$ defined in (2.3). From this viewpoint, the unadjusted estimator $\widehat{\tau}_{1}^{\text{UNADJ}}$ also belongs to the same class of balancing weighting estimators as $\widehat{\tau}_{1}^{h}$. Therefore, we can construct the influence function of $\widehat{\tau}_{1}^{\text{UNADJ}}$ (denoted as $\phi^{\text{UNADJ}}(\mathcal{O})$ — see Web Appendix B for details) in the format of (3.7) but without the second term on the right hand side. Typically, for the unadjusted estimator $\widehat{\tau}_{1}^{\text{UNADJ}}$, $e(X)=\pi$ and $\tau_{1}^{h}=\tau_{1}$; then we can show $\phi^{\text{UNADJ}}(\mathcal{O})=\chi^{h}(\mathcal{O})$ and ${\text{var}}\left(\widehat{\tau}_{1}^{\text{UNADJ}}\right)={\text{var}}\left(\chi^{h}(\mathcal{O})\right)$. Theorem 1(c) shows that the influence function of a balancing weighting estimator can be expressed as the influence function of the unadjusted estimator minus its projection onto the nuisance tangent space, which in fact underlies the asymptotic variance reduction achieved through covariate adjustment by weighting in randomized trials.

We note that Theorem 1 summarizes the large-sample properties of the IPW and OW estimators for win estimands in randomized trials, which explains the asymptotic variance reduction of weighting estimators theoretically, extending the findings demonstrated in Shen et al. (2014), Williamson et al. (2014) and Zeng et al. (2021). However, the weighting estimators derived in this article are based on U-statistics and our win estimand is more complex than the usual average treatment effect.

Remark** 2**

The large-sample properties of the OW estimator derived in this article do not depend on equal randomization ($\pi=0.5$), which is different from that in Zeng et al. (2021). The reason is that the OW estimator is constructed by a symmetric kernel function based on two weights, such that the expectation of the kernel function of OW is equal to that of IPW under large sample sizes (i.e., $\mathbb{E}(q^{\text{IPW}}(\mathcal{O}_{i},\mathcal{O}_{j};\widehat{\beta}))=\mathbb{E}(q^{\text{OW}}(\mathcal{O}_{i},\mathcal{O}_{j};\widehat{\beta}))$). Therefore, equal randomization is not necessary for the variance reduction property under either OW or IPW.

3.4 Closed-form variance estimation

By Theorem 1(a), ${\text{var}}(\widehat{\tau}_{1}^{h})$ can be consistently estimated as

[TABLE]

where $\widehat{q}_{1}^{h}(\mathcal{O}_{i};\widehat{\beta})=(n-1)^{-1}\sum_{j\neq i}^{n}q^{h}\left(\mathcal{O}_{i},\mathcal{O}_{j};\widehat{\beta}\right)$ and $A_{1}^{h}={n\choose 2}^{-1}\sum_{i=1}^{n}\sum_{j>i}^{n}\frac{\partial}{\partial\beta}q^{h}\left(\mathcal{O}_{i},\mathcal{O}_{j};\beta\right)|_{\beta=\widehat{\beta}}$, $\widetilde{l}_{\beta}(Z,X;\widehat{\beta})$ is the influence function for $\widehat{\beta}$ and $a^{\otimes 2}=aa^{\top}$ for any vector $a$. Therefore, for $\widehat{\tau}_{1}^{\text{IPW}}$ (i.e., $\widehat{\tau}_{1}^{h}$ with $\omega(X_{i},X_{j})=1/\widehat{e}(X_{i})[1-\widehat{e}(X_{j})]$) and $\widehat{\tau}_{1}^{\text{OW}}$ (i.e., $\widehat{\tau}_{1}^{h}$ with $\omega(X_{i},X_{j})=\widehat{e}(X_{j})[1-\widehat{e}(X_{i})]$),we can obtain their consistent variance estimators based on (3.8), and their explicit expressions are provided in Web Appendix C.

In fact, according to the result in Theorem 1(c) and (3.8), the variance of $\widehat{\tau}_{1}^{\text{UNADJ}}$ can also be consistently estimated as the second moment estimator for an estimated influence function:

[TABLE]

where $\widehat{q}_{1}^{\text{UNADJ}}(\mathcal{O}_{i})=(n-1)^{-1}\sum_{j\neq i}^{n}q^{\text{UNADJ}}(\mathcal{O}_{i},\mathcal{O}_{j})$ and $q^{\text{UNADJ}}(\mathcal{O}_{i},\mathcal{O}_{j})=\frac{1}{2}\bigr[Z_{i}(1-Z_{j})\mathbb{I}(Y_{i}>Y_{j})+Z_{j}(1-Z_{i})\mathbb{I}(Y_{i}<Y_{j})\bigr]/\bigr[\frac{1}{n(n-1)}\sum_{i=1}^{n}\sum_{j=1}^{n}Z_{i}(1-Z_{j})\bigr]$. In Section 5, we will demonstrate that the variance estimators (3.9) and (2.5) perform similarly and are both valid when sample size is moderate, but that (3.9) generally has better small sample performance than (2.5).

Finally, since we have established the asymptotic results for a general kernel function $g$, the results would be applicable for both the win and loss estimators and can inform the estimation of WR and WD through application of the delta method. Specifically, for the loss probability estimator

[TABLE]

we can obtain the corresponding variance estimators by simply changing the kernel function $g$.

Finally, denote $\phi_{1,i}^{h}=2\left\{\widehat{q}_{1}^{h}(\mathcal{O}_{i};\widehat{\beta})-\widehat{\tau}_{1}^{h}\right\}+\widetilde{l}^{\top}_{\beta}(Z_{i},X_{i};\widehat{\beta})A_{1}^{h}$ and $\phi_{-1,i}^{h}=2\left\{\widehat{q}_{-1}^{h}(\mathcal{O}_{i};\widehat{\beta})-\widehat{\tau}_{-1}^{h}\right\}+\widetilde{l}^{\top}_{\beta}(Z_{i},X_{i};\widehat{\beta})A_{-1}^{h}$ as the influence function of patient $i$ for the probability of win and probability of loss, respectively, where $\widehat{q}_{-1}^{h}(\mathcal{O}_{i};\widehat{\beta})$, $\widehat{\tau}_{-1}^{h}$ and $A_{-1}^{h}$ are the corresponding estimators for the probability of loss — see (3.8) for similar definitions. Then, the asymptotic covariance between $\widehat{\tau}_{1}^{h}$ and $\widehat{\tau}_{-1}^{h}$ can be computed as ${\text{cov}}\left(\widehat{\tau}_{1}^{h},\widehat{\tau}_{-1}^{h}\right)=n^{-1}\sum_{i=1}^{n}\phi_{1,i}^{h}\phi_{-1,i}^{h}$. By applying formulas (2.6) and (2.7), we can then assemble the asymptotic variance and covariance expressions to obtain final variance estimators for WR ($\theta$) and WD ($\eta$) based on IPW and OW.

4 Covariate-adjusted win statistics via augmented weighting

4.1 Formulating an augmented IPW estimator with outcome regression

Both IPW and OW estimators only use a working propensity score model to adjust for chance imbalance. Hence, intuitively, potentially more efficient estimators could be constructed by leveraging the associations between baseline covariates and the potential outcomes through an additional outcome regression model, that is, the augmented weighting methods. For observational studies, Mao (2018) proposed locally efficient and doubly robust estimators for general causal estimands based on contrast functions, formulated using Horvitz-Thompson weighting. Mao et al. (2019) extended this framework to study augmented weighting estimators, including augmented IPW and OW variants for weighted average treatment effects, and derived closed-form variance expressions. More recently, Buyse et al. (2025) discussed an augmented IPW estimator for generalized pairwise comparison estimands, using Horvitz–Thompson weights. In what follows, we first review the formulation of the augmented IPW estimator discussed in Mao (2018). Building on this estimator, we then extend the approach to construct stabilized augmented weighting estimators for the win estimand allowing for the general class of balancing weights.

Let $\nu_{z}(\cdot|X;\alpha)$ denote the conditional laws (i.e., conditional density for continuous random variable or conditional probability for discrete random variable) of $Y(Z)$ given $X$, indexed by a Euclidean parameter $\alpha$, and let $\mathcal{M}_{\mu}$ denote the intersection of these two marginal models. For notational simplicity, we assume both $\mathcal{M}_{\mu}$ and $\mathcal{M}_{e}$ use the same set of covariates, though in practice these may differ. Given two independent individuals with $X_{i}$ and $X_{j}$, an outcome regression component can be defined as

[TABLE]

After finding the maximum likelihood estimator $\widehat{\alpha}$ for $\alpha$ under $\mathcal{M}_{\mu}$ based on $(Y_{i},Z_{i},X_{i})(i=1,\dots,n)$, and obtaining $\widehat{\mu}(X_{i},X_{j})=\mu(X_{i},X_{j};\widehat{\alpha})$, the augmented IPW estimator for $\tau^{g}$ in Mao (2018) is given by

[TABLE]

where

[TABLE]

Alternatively, we can rewrite the above augmented weighting estimator as follows:

[TABLE]

Specifically, for an ordinal outcome $Y$, let $p_{l}^{z}(X;\alpha)$ denote the model-based conditional probability that $Y(Z)=l$ $(l=1,2,\dots,L)$ given baseline variables $X$ and parameter $\alpha$, then $p_{l}^{z}(X;\alpha)=\nu_{z}(Y(Z)=l|X;\alpha)$. Assume that individual $i$ belongs to the treatment group ($Z=1$) and individual $j$ to the control group ($Z=0$). Then, Mao (2018) has shown that

[TABLE]

In our simulation study, we primarily focus on ordinal logistic regression to model $\mu(\cdot,\cdot;\alpha)$ and to obtain fitted estimates using the polr function from the R package MASS (Venables and Ripley, 2002). In the polr function, the ordinal logistic regression model for an ordinal outcome with $L$ categories is parameterized as

[TABLE]

where $l=1,2,\dots,L-1$, $\xi=(\xi_{1},\dots,\xi_{p})^{\top}$ and $\alpha=(\alpha_{10},\dots,\alpha_{L-1,0},\xi^{\top})^{\top}$. Then, we have

[TABLE]

However, the augmented IPW estimator is fairly general and in theory can be used with any working outcome models, including more flexible ordinal regression models that permit the odds ratios for covariates to possibly change with cut-point.

4.2 Covariate adjustment via augmented inverse probability weighting and overlap weighting

It is worth noting that when using the augmented IPW estimator $\widehat{\tau}_{1}^{\text{MAO}}$, the pairwise comparison concerning two individuals from the same treatment group (i.e., $Z_{i}=Z_{j}=1$ or $Z_{i}=Z_{j}=0$) simplifies to $r^{\text{MAO}}(\mathcal{O}_{i},\mathcal{O}_{j};\widehat{\beta},\widehat{\alpha})=\bigr[\widehat{\mu}(X_{i},X_{j})+\widehat{\mu}(X_{j},X_{i})\bigr]/2$. Therefore, the estimator $\widehat{\tau}^{\text{MAO}}_{1}$ can be reexpressed as

[TABLE]

For an ordinal outcome, the kernel function used to define the win probability $\tau_{1}$ in (2.1) is given by $g(Y_{i},Y_{j})=\mathbb{I}(Y_{i}>Y_{j})$. Since the estimator $\tilde{\tau}_{1}^{\text{MAO}}$ is not a Hájek-type estimator, we incorporate a tilting function $\widehat{h}(X_{i},X_{j})$ and stabilized weights to propose the following augmented weighting (AW) estimator:

[TABLE]

where $\omega(X_{i},X_{j})=\widehat{h}(X_{i},X_{j})/\{\widehat{e}(X_{i})(1-\widehat{e}(X_{j}))\}$. In the case where IPW is used, the tilting function is $\widehat{h}(X_{i},X_{j})=1$ and the weight function is defined as $\omega(X_{i},X_{j})=\{\widehat{e}(X_{i})(1-\widehat{e}(X_{j}))\}^{-1}$. This leads to the following augmented inverse probability weighting (AIPW) estimator:

[TABLE]

Alternatively, by setting the tilting function to $\widehat{h}(X_{i},X_{j})=\widehat{e}(X_{i})\widehat{e}(X_{j})(1-\widehat{e}(X_{i}))(1-\widehat{e}(X_{j}))$, we have $\omega(X_{i},X_{j})=\{1-\widehat{e}(X_{i})\}\widehat{e}(X_{j})$ and the following augmented overlap weighting (AOW) estimator is defined as:

[TABLE]

4.3 Asymptotic analysis of the AIPW and AOW estimators for win estimands

To derive the variance of $\widehat{\tau}_{1}^{\text{AW}}$ as defined in (4.1), we first express it in the form of a U-statistic:

[TABLE]

where

[TABLE]

By the uniform law of large numbers for U-processes (Arcones and Giné, 1993; Van Der Vaart and Wellner, 1996) and the consistency of both $\widehat{\beta}$ and $\widehat{\alpha}$, $\widehat{\tau}_{1}^{\text{AW}}$ is consistent when either $\mathcal{M}_{e}$ or $\mathcal{M}_{\mu}$ is correctly specified. Theorem 2 below provides the asymptotic variance of the augmented weighting estimator and establishes the local efficiency and robustness of $\widehat{\tau}_{1}^{\text{AW}}$.

Theory** 2**

*(Asymptotic results of augmented weighting estimators)

(a) The asymptotic variance of $n^{1/2}(\widehat{\tau}_{1}^{\text{AW}}-\tau_{1}^{h})$ is ${\text{var}}\left(\psi^{h}(\mathcal{O})\right)$, where*

[TABLE]

is the influence function of $\widehat{\tau}_{1}^{\text{AW}}$ with

[TABLE]

*and $\widetilde{l}_{\alpha}\left(\mathcal{O};\alpha\right)$ is the influence function for $\alpha$.

(b) Under randomization, both the AIPW and AOW estimators, $\widehat{\tau}_{1}^{\text{AIPW}}$ and $\widehat{\tau}_{1}^{\text{AOW}}$ are robust in the sense that they are always consistent to $\tau_{1}^{h}$ regardless of whether the ordinal outcome regression model $\mathcal{M}_{\mu}$ is correctly specified. When the ordinal outcome regression model $\mathcal{M}_{\mu}$ is correctly specified, both estimators are locally efficient.***

The consistent variance estimators for $\widehat{\tau}_{1}^{\text{AW}}$ can be constructed similarly to those for $\widehat{\tau}_{1}^{h}$ in (3.8). According to lemma S1 of Mao (2018), we have $2[E\{r^{h}(\mathcal{O}_{i},\mathcal{O}_{j};\beta,\alpha)|\mathcal{O}_{i}\}-\tau_{1}^{h}]=\psi^{h}(\mathcal{O}_{i})$, then by Hoeffding’s decomposition theorem (Hoeffding, 1948), the remainder of the proof for Theorem 2(b) follows the same steps as for the average treatment effect discussed in Chapter 13 of Tsiatis (2006). Typically, when $\mathcal{M}_{\mu}$ is correct, $\widehat{\alpha}$ has no effect on the asymptotic distribution of $\widehat{\tau}_{1}^{\text{AW}}$. In other words, $\widetilde{l}_{\alpha}\left(\mathcal{O};\alpha\right)$ has no contribution to $\psi^{h}(\mathcal{O})$ and $\psi^{h}(\mathcal{O})$ can be rewritten as

[TABLE]

Note that for $\tau_{1}^{\text{UNADJ}}$, $\pi=e(X)$, so $2\{r_{1}^{h}(\mathcal{O};\beta,\alpha)-\tau_{1}^{h}\}=\chi^{h}(\mathcal{O})$ (i.e., the influence function of $\tau_{1}^{\text{UNADJ}}$). By using the same technique in Theorem 1(b), we can also show that $\psi^{h}(\mathcal{O})$ is equal to $\chi^{h}(\mathcal{O})$ minus its projection on the tangent space of $\mathcal{M}_{e}$. Therefore, the asymptotic variance reduction by augmented weighting estimators is guaranteed when the outcome model is correctly specified, a finding that is consistent to that in Vermeulen et al. (2015) for studying covariate-adjusted Mann-Whitney U tests.

Finally, based on Theorem 2, we can also obtain the closed-form variance estimators of $\widehat{\tau}_{1}^{\text{AIPW}}$ and $\widehat{\tau}_{1}^{\text{AOW}}$; the technical details and full expressions are provided in Web Appendix D. Similarly, we can apply the procedure described in the penultimate paragraph of Section 3.4 to derive the covariance between the win probability $\widehat{\tau}_{1}^{\text{AW}}$ and loss probability $\widehat{\tau}_{-1}^{\text{AW}}$. This allows us to estimate the variance of $\widehat{\theta}$ using formula (2.6), and the variance of $\widehat{\eta}$ using formula (2.7), for both the AIPW and AOW estimators. To assist practical applications, we implement the unadjusted and all covariate-adjusted point and variance estimators in an R package winPSW.

5 Simulation studies

To investigate the finite-sample performance of unadjusted, weighting and augmented weighting estimators, we conduct a series of simulation studies under several data generation processes. Our overarching goal is to empirically (1) examine the accuracy of the point estimators for WR and WD with ordinal outcomes; (2) validate the variance estimators of the unadjusted, weighting and augmented weighting methods developed in (2.5) or (3.9), (A.8), (A.9) in Web Appendix C, (A.13) and (A.14) in Web Appendix D, respectively; (3) explore whether the weighting and augmented weighting estimators can lead to finite-sample efficiency gains compared with the unadjusted estimator. Through these simulation studies, we also seek to confirm the model robustness property of augmented weighting estimators under a mis-specified working ordinal outcome model.

5.1 Simulation design

We generate $p=6$ baseline covariates (three continuous and three binary covariates): $X_{1}\sim N(1,0.3^{2})$, $X_{2}\sim N(0.9,0.4^{2})$, $X_{3}\sim N(0.8,0.5^{2})$, $X_{4}\sim\text{Bernoulli}(0.75)$, $X_{5}\sim\text{Bernoulli}(0.5)$, and $X_{6}\sim\text{Bernoulli}(0.25)$. The treatment indicator is generated from $Z\sim\text{Bernoulli}(\pi)$ by fixing the randomization probability $\pi$. For illustrative purposes, we assume the ordinal outcome has three levels and that a higher level indicates a more favorable outcome. To examine the model robustness feature of our covariate-adjusted estimators, we consider two types of potential outcomes data generating processes: a quadratic association cumulative logit model, and a covariate interaction cumulative logit model. Specifically, based on baseline covariates $X_{i}=(X_{i1},X_{i2},\dots,X_{i6})^{\top}$, we generate the potential outcomes from the following quadratic model: for $z=0,1$,

[TABLE]

where $\alpha_{10}=1$, $\alpha_{20}=0.05$, $\delta=1$, $\gamma=(\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4},\gamma_{5},\gamma_{6})^{\top}=(1,-1,1,-1,1,-1)^{\top}\times(0.5+0.5z)$ and $\gamma_{quad}=(\gamma_{quad,1},\gamma_{quad,2},\gamma_{quad,3},\gamma_{quad,4},\gamma_{quad,5},\gamma_{quad,6})^{\top}=2\gamma$. Therefore, the parameter $\alpha$ in outcome $\mu(X_{i},X_{j};\alpha)$ can be characterized as $\alpha=(\alpha_{10},\alpha_{20},\gamma^{\top},\gamma^{\top}_{quad},\delta)^{\top}$. In addition, we consider the following interaction model to generate the potential outcomes: for $z=0,1$,

[TABLE]

where $\gamma_{int,jk}=\mathbb{I}(\gamma_{j}\times\gamma_{k}>0)\times(\gamma_{j}+\gamma_{k})+\mathbb{I}(\gamma_{j}\times\gamma_{k}<0)\times 0.25\times(\gamma_{j}\times\gamma_{k})$ and $\delta=2$. Values of other parameters $\alpha_{10}$, $\alpha_{20}$ and $\gamma$ are the same as those in quadratic model (5.1). The observed outcome is set to be $Y_{i}=Y_{i}(Z_{i})=Z_{i}Y_{i}(1)+(1-Z_{i})Y_{i}(0)$. For the randomization probability $\pi$, two values are considered: $\pi=0.5$ and $\pi=0.7$, representing a balanced assignment and an unbalanced assignment.

Under each simulation scenario, true values of WR $\theta$ and WD $\eta$ are determined using Monte Carlo approximation based on the simulated potential outcomes. We compare the unadjusted estimator $\widehat{\tau}^{\text{UNADJ}}$ (UNADJ), the IPW estimator $\widehat{\tau}^{\text{IPW}}$ (IPW), the OW estimator $\widehat{\tau}^{\text{OW}}$ (OW), the augmented IPW estimator $\widehat{\tau}^{\text{AIPW}}$ (AIPW) and the augmented OW estimator $\widehat{\tau}^{\text{AOW}}$ (AOW) for WR and WD. For weighting estimators and augmented weighting estimators, the propensity score is estimated by logistic regression including six baseline covariates as linear terms. For the augmented weighting estimators, we keep the same propensity score model specification, and the working outcome regression model is specified as an ordinal logistic regression including the six baseline covariates as linear terms plus quadratic terms or plus interaction terms.

When evaluating the accuracy of the point estimator, we generate the observed data $\mathcal{O}_{i}=(Y_{i},Z_{i},X_{i})$ with three different sample sizes $n\in\{200,300,400\}$ under both balanced and unbalanced designs based on the quadratic association model (5.1) and the covariate interaction model (5.2), and calculate the relative bias (bias divided by the true estimand value times 100). To explore the model robustness property of the augmented weighting estimators, we also explore the results when we omit quadratic terms $X^{2}_{i,j}$ and interaction terms $X_{i,j}X_{i,k}$ ($j\neq k$) when specifying the ordinal outcome logistic regression. We refer to this scenario as “model mis-specification” and corresponding augmented IPW and augmented OW as AIPW-mis and AOW-mis, respectively.

Beyond the accuracy of the point estimators, we also investigate the accuracy of the proposed variance estimators and the empirical coverage of the 95% confidence intervals (based on normal approximation). Specifically, for $k=1$ or $k=-1$, the confidence intervals of $\widehat{\tau}_{k}^{\text{UNADJ}}$, $\widehat{\tau}_{k}^{\text{IPW}}$, $\widehat{\tau}_{k}^{\text{OW}}$, $\widehat{\tau}_{k}^{\text{AIPW}}$, $\widehat{\tau}_{k}^{\text{AOW}}$ are constructed based on (3.9), (A.8), (A.9), (A.13) and (A.14) in Web Appendix, respectively. Then, based on (2.6) and (2.7), we can obtain confidence intervals of $\widehat{\theta}$ and $\widehat{\eta}$. To quantify the accuracy of the estimated variance of each estimator, we report the ratio between the average variance estimates and the Monte Carlo variance (considered as the truth) of each estimator, with a ratio close to 1 indicating adequate performance. In addition, we also compare the variance estimator of $\widehat{\tau}_{1}^{\text{UNADJ}}$ based on (2.5) (the original Bebu and Lachin (2016) approach), the bias-reduced version provided in Web Appendix A3, and (3.9) (our alternative variance estimation approach). We report relative bias (%), Monte Carlo standard deviation (MCSD), average estimated standard error (AESE) based on (2.5) and (3.9), and corresponding coverage probability (CP) of the 95% confidence intervals for the WR and WD.

To illustrate the finite-sample efficiency gains under different sample sizes, we designed a separate set of simulations based on both the quadratic model (5.1) and the interaction model (5.2), where the total sample sizes $n$ vary from 60 to 600 to mimic small and large sample size scenarios. We report the efficiency results for estimating WR and WD based on the Monte Carlo variance across all simulation scenarios. The relative efficiency of an estimator $\widehat{\theta}$ (or $\widehat{\eta}$) is defined as ${{\text{var}}(\widehat{\theta}^{\text{UNADJ}})}/{{\text{var}}(\widehat{\theta})}$ (or ${{\text{var}}(\widehat{\eta}^{\text{UNADJ}})}/{{\text{var}}(\widehat{\eta})}$), i.e., the ratio between the Monte Carlo variance of the unadjusted estimator and that of the estimator under consideration. A higher relative efficiency indicates that the estimator is more efficient than the unadjusted estimator, and hence covariate adjustment is beneficial. We also report the relative efficiency comparisons between the augmented weighting estimators and the unadjusted estimator when the working outcome regression model is mis-specified (i.e., the quadratic terms and interaction terms are omitted when fitting the outcome $Y$ by an ordinal logistic regression). For each scenario, we simulate 1000 data replicates and the results are summarized based on these replicates.

5.2 Simulation results

5.2.1 Point estimation

[Table 1 about here.]

The relative bias of UNADJ, IPW, OW, AIPW, AOW, AIPW-mis and AOW-mis under each scenario are summarized in Table 1. For the WR, we notice slight relative bias from all estimators even when the sample size is as small as 200. When the sample size increases to 400, all relative biases for estimating WR from weighting and augmented weighting estimators are negligible (often less than 1%). Furthermore, under each scenario, the bias of the WR from the unadjusted estimator is usually the largest, suggesting that covariate adjustment may reduce finite-sample bias. For estimating the WD, the relative biases from all estimators remain negligible and the bias from the unadjusted estimator is again slightly larger than that from the covariate-adjusted estimators. However, under an unbalanced design and when sample size is as small as 200, the bias from the unadjusted estimator can be slightly smaller than that from the covariate-adjusted estimators. Comparing biases of AIPW-mis and AOW-mis with those of AIPW and AOW, we find that the relative bias of AIPW and AOW is not substantially different regardless of whether the outcome model is correctly specified. This highlights the outcome model robustness property of the augmented weighting estimators.

5.2.2 Variance and interval estimation

[Table 2 about here.]

Table 2 reports the results on the accuracy of the proposed variance estimators and the empirical coverage rates of the associated 95% confidence interval estimators for WR and WD under both balanced and unbalanced designs. Overall, we find that the average estimated variance to Monte Carlo variance ratio is generally close to 1 under all scenarios, especially when the sample size reaches $n=400$. This pattern clearly demonstrates that all of the proposed variance estimators are generally accurate in capturing the true variability of the covariate-adjusted estimators. Noticeably, even when sample size reaches 400, the ratio between the average estimated variance to Monte Carlo variance from both AIPW and AOW under the covariate interaction data generating model is generally less than 0.9; this is especially the case for AIPW under the unbalanced design, where this ratio drops below 0.8. Upon closer inspection, we find that when fitting both AIPW and AOW under the correct outcome model, there are a handful of baseline covariates (6 baseline covariates plus 15 interaction terms) in the fitted outcome model. Hence, the increasing number of coefficients that need to be estimated compromises the finite-sample convergence of the variance estimator to the truth. When the misspecified outcome model only includes 6 baseline covariates (AIPW-mis and AOW-mis), the ratio between the average estimated variance to Monte Carlo variance are closer to 1. This result shows that for augmented weighting estimators, if the sample size is not large, including too many covariates in the outcome model may lead to underestimation of variance. Finally, as the sample size increases, coverage rates generally improve across all estimators, reflecting better estimation of the true variance. Particularly, our proposed variance estimators accurately capture the true variability of AIPW-mis and AOW-mis even under unbalanced designs.

We also explore the performance of the proposed variance estimator against the original Bebu and Lachin (2016) variance estimator and the bias-reduced variance estimator shown in Web Appendix A3 (suggested by Ozenne et al. (2021)) for the UNADJ estimator. Web Table S1 shows that the average estimated standard error calculated by (3.9) (AESE1), (2.5) (AESE2) and the method developed in Web Appendix A3 (AESE3) are close to the Monte Carlo standard deviation (MCSD) for all estimators when sample size $n\geq 200$, and they are associated with nominal coverage, especially when the sample size is $n=400$. Furthermore, AESE1 and AESE3 are generally close, which is expected since these two corresponding variance estimators share the same underlying mechanism (projection theorem of U-statistics). However, with a small sample size (e.g., $n\leq 50$), coverage probabilities associated with AESE1 and AESE3 are closer to nominal than that associated with AESE2, especially under an unbalanced design. This result confirms that, if no covariate adjustment is performed, all three variance estimators are reliable when sample size is relatively moderate (e.g., $n\geq 200$), but only (3.9) and bias-reduced variance estimator in Web Appendix A3 are recommended for smaller sample sizes.

5.2.3 Relative efficiency

Figures 1-2 present the relative efficiency of IPW, OW, AIPW and AOW (AIPW-mis and AOW-mis) for WR and WD across sample sizes ranging from 60 to 600 when the data generating process follows the quadratic association model (Figure 1) and the covariate interaction model (Figure 2), respectively. In each figure, we consider four scenarios: (A) balanced design with a correctly specified outcome regression model; (B) balanced design with a mis-specified outcome regression model; (C) unbalanced design with a correctly specified outcome regression model; (D) unbalanced design with a mis-specified outcome regression model.

[Figure 1 about here.]

[Figure 2 about here.]

First, the propensity score estimators are notably more efficient than the unadjusted estimator, which confirms that adjusting for prognostic covariates can significantly improve the efficiency for estimating WR and WD, sometimes by as much as 2.5-fold. Furthermore, the augmented weighting estimators coupled with an additional outcome regression model (correctly or incorrectly specified) generally show additional efficiency gains across all scenarios. The only exceptions are (i) when the sample size is relatively small ($n\leq 120$) under an unbalanced design (row C in Figure 1), and (ii) when the data generating process follows the covariate interaction model under an unbalanced design (row C in Figure 2). The likely reason for the latter case is that, in the true covariate interaction model (5.2), there are 15 interaction terms in the true outcome model, posing stronger sample size requirements for the correctly specified AIPW and AOW estimators to gain further efficiency over the simpler propensity score weighting estimators. Except for that challenging scenario, when mis-specification of the outcome regression is present, the relative efficiency of AIPW and AOW is not much different than the case where the outcome regression is correctly specified, especially when the sample size is relative large, say $n\geq 420$ (see panel A vs panel B, panel C vs panel D in Figures 1-2). Finally, when mis-specification occurs in the scenario of an unbalanced design and the covariate interaction data generating model, the augmented weighting estimators still lead to higher efficiency over the propensity score weighting estimators (see panel D in Figure 2) over the entire range of sample sizes. This suggests that, in relatively small RCTs with unbalanced assignment, careful consideration should be given to outcome model specification when using AIPW or AOW estimators, particularly regarding the number of covariates to adjust for, due to the tradeoff between model complexity and limited sample size.

These simulation results also shed light on the efficiency comparison between using IPW and OW for covariate adjustment when estimating WR and WD. Under a balanced design, we find that IPW and OW, as well as AIPW and AOW often perform similarly in terms of efficiency, especially as the sample size increases. The use of overlap weights can have a slight efficiency advantage when the sample size is small, say, no more than $120$. The efficiency difference between the use of different weights (IPW versus OW, AIPW versus OW) is more obvious under an unbalanced design. Under an unbalanced design which inherently offers fewer pairwise comparisons, OW-based estimators often demonstrate superior efficiency over their IPW-based counterparts, especially when sample sizes are smaller. Furthermore, the efficiency advantages of AOW over AIPW and OW over IPW are more pronounced when estimating the WD across all scenarios within unbalanced designs and when estimating the WR in scenarios under the quadratic association data generating model within relatively small sample size (e.g., $n\leq 120$) unbalanced design. Thus, the relative efficiency of OW and AOW estimators is most beneficial in smaller studies or those with unbalanced randomization.

5.3 Additional simulations for the testing procedures

Theorem 1 and 2 also provide a principled basis for deriving testing procedures for $\Delta\in\{\log(\theta),\eta\}$. Consider the Wald-type test based on the studentized statistic $Z_{*}=\widehat{\Delta}/\widehat{\sigma}_{*}$, where $\widehat{\sigma}_{*}$ consistently estimates the true standard deviation $\sigma_{*}$. Then under the null hypothesis $H_{0}:\Delta=0$, we have $Z_{*}\stackrel{{\scriptstyle d}}{{\sim}}N(0,1)$. Note that for both model (5.1) and model (5.2), if we let the regression coefficient $\delta=0$, $Y_{i}(1)$ and $Y_{j}(0)$ are generated from the exact same distribution, allowing us to explore type I error of UNADJ, IPW, OW, AIPW, AOW, AIPW-mis and AOW-mis. The empirical type I error rate results under Wald $z$ test are reported in Table 3. We observe that under the null, the empirical type I error rates from UNADJ, IPW, OW, AOW, AIPW-mis and AOW-mis generally remain close to the nominal 5% level, but the empirical type I error rate from AIPW can be inflated and this is more pronounced for a small sample and an unbalanced design. This result is consistent with the simulation findings in Section LABEL:mc_var and is likely due to the complexity of the true outcome model. That is, when the outcome regression includes many interaction terms relative to the available sample size, estimation of the augmentation component and its associated variance can become less stable, leading to mild finite-sample inflation of the Wald test.

[Table 3 about here.]

Further, we conducted additional simulations to evaluate whether results in Section LABEL:mc_bias and LABEL:mc_var will be affected by strong correlation among covariates or not. Detailed data generation when covariates are strongly correlated is presented in Web Appendix E, while corresponding point estimates, variance and interval estimates, and empirical type I error rate are shown in Tables S2, S3 and S4, respectively. Since results in Table S2 and Table 1, Table S3 and Table 2, as well as Table S4 and Table 3, are similar and comparable, it is concluded that our findings based on Section 5.2 are not affected by strong correlation among covariates.

6 An Illustrative Application to the ORCHID clinical trial

To illustrate our proposed methods, we reanalyze a completed RCT — the Outcomes Related to COVID-19 Treated With Hydroxychloroquine Among Inpatients With Symptomatic Disease (ORCHID). This multi-center, double-blind, placebo-controlled trial was carried out under the PETAL Clinical Trials Network between April and June 2020, and aims to test the hypothesis that, compared with placebo, hydroxychloroquine improves clinical outcomes for adults hospitalized with COVID-19 (Self et al., 2020). The study enrolled 479 adults hospitalized with confirmed SARS-CoV-2 infection and respiratory symptoms, and randomly assigned them to receive either hydroxychloroquine or placebo. The primary outcome was the patient’s clinical status at 14 days, assessed via a 7-point ordinal scale: (1) Not hospitalized, performing normal activities; (2) Not hospitalized, but unable to perform normal activities; (3) Hospitalized, not requiring supplemental oxygen; (4) Hospitalized, requiring supplemental oxygen via nasal cannula; (5) Hospitalized, requiring respiratory support via non-rebreather mask, or high-flow nasal cannula oxygen therapy, or noninvasive mechanical ventilation; (6) Hospitalized, requiring extracorporeal membrane oxygen, invasive mechanical ventilation, or both; (7) Death; all of which are based on the World Health Organization-recommended 7-point COVID Outcomes Scale (Cummings et al., 2020). The trial also included 12 secondary outcomes, including 7-day and 28-day mortality outcomes. Despite early optimism based on preclinical antiviral and immunomodulatory properties of hydroxychloroquine, the primary analysis based on a multivariable proportional odds model obtained an adjusted odds ratio (aOR, a value greater than 1.0 indicating more favorable outcomes with hydroxychloroquine than placebo) of 1.02, with 95% confidence interval [0.73, 1.42]. This result means that, compared with placebo, treatment with hydroxychloroquine did not significantly improve clinical status at day 14. We note that this estimand is defined through the regression model rather than as a marginal contrast of potential outcomes, and that it relies on the proportional odds assumption; in contrast, the win estimands we will consider below are marginal, model-free summaries based on pairwise comparisons of potential outcomes, providing a complementary perspective on the same data.

In this data example, we leverage our win statistics framework to analyze the ordinal clinical status outcomes measured at 7, 14, and 28 days after randomization, with the 14-day outcome corresponding to the primary outcome in the original ORCHID trial analysis. We adjust for a comprehensive set of potential prognostic factors, including body mass index, ICU status at baseline, standardized age, race, sex, baseline oxygen saturation, respiratory rate, and the Glasgow Coma Scale component of the SOFA score. These covariates are adjusted for as linear terms in the logistic propensity score model for IPW and OW estimators. Moreover, for augmented weighting estimators, we also adjust for the same set of covariates as linear terms in the ordinal logistic outcome model $\mu(\cdot,\cdot;\alpha)$. To quantify the treatment effect of hydroxychloroquine compared to placebo, we estimate the WR and WD with the five different methods (UNADJ, IPW, OW, AIPW and AOW). The point estimate, standard error (SE) calculated by corresponding closed-form variance estimation ((3.9), (A.8), (A.9), (A.13) and (A.14) in Web Appendix, respectively) as well as the 95% confidence intervals are summarized in Table 4. In addition, to measure the finite-sample variance reduction from weighting and augmented weighting estimators (compared to the unadjusted estimator) through covariate adjustment, we compute the proportion of variance reduction compared to the unadjusted estimator (PVR), defined as $\text{PVR}={(\text{SE}^{2}_{\text{UNADJ}}-\text{SE}^{2}_{\text{Estimator}})}/{\text{SE}^{2}_{\text{UNADJ}}}\times 100\%$, where $\text{Estimator}\in\{\text{IPW, OW, AIPW, AOW}\}$.

[Table 4 about here.]

According to Table 4, we observe that across all three time points (Day 7, Day 14, and Day 28), the estimated win ratios for hydroxychloroquine versus placebo from the five proposed methods carry magnitudes larger than 1, with 95% confidence intervals consistently including the null, indicating no statistically significant treatment effect on the pairwise-comparison scale. The direction and statistical conclusions when the win difference is reported are the same. Both estimated WR and WD conclusions are consistent with the original analysis in Self et al. (2020). To interpret the results, we consider the AOW estimates for the outcome at Day 7 as an example. Under this approach, the WR comparing hydroxychloroquine to placebo was 1.066 (95% CI: 0.855–1.278), indicating that the probability of a more favorable outcome for a randomly selected hydroxychloroquine-treated patient compared to a placebo patient was 6.6% higher than the reverse. Similarly, the WD comparing hydroxychloroquine to placebo was 0.026 (95% CI: -0.055–0.108), indicating that the probability of a more favorable outcome for a randomly selected hydroxychloroquine-treated patient exceeded that for a placebo patient by 2.6 percentage points. Finally, comparing across methods, the point estimates are generally similar regardless of whether baseline covariate adjustment is used. However, the standard error estimates after covariate adjustment are typically smaller. For example, when estimating the WR for the Day 14 outcome, covariate-adjusted weighting estimators (IPW and OW) and augmented weighting estimators (AIPW and AOW) achieve over 30% variance reduction (34.4%, 33.3%, 34.1% and 33.5%, respectively) compared with unadjusted estimators. Overall, we find that the amount of variance reduction is relatively consistent across the weighting and augmented weighting estimators, suggesting that including an additional ordinal outcome regression beyond weighting, in this example, does not further reduce the variance of the WR and WD estimators.

7 Discussion

In this article, we have introduced covariate-adjustment methods to estimate the win ratio and win difference. We first developed propensity score weighting estimators — in particular, IPW and OW — for win statistics. These estimators extend the prior work of Zeng et al. (2021) from the conventional average treatment effect estimand to a wider class of pairwise comparison estimands, for which win ratios and win differences are important special cases. Second, inspired by the form of the augmented weighting estimator in Mao et al. (2019) for analyzing observational data, we further introduced a class of augmented weighting estimators for win estimands that leveraged an additional outcome regression (in addition to weighting) to potentially further improve efficiency. Both weighting and augmented weighting estimators (IPW or OW) were constructed using stabilized weights such that they are Hájek-type estimators, different from the unstabilized IPW and augmented inverse probability weighting estimators previously proposed for observational studies (Mao, 2018). We applied the U-statistics theory to derive the asymptotic variance of all unadjusted and covariate-adjusted estimators, and obtained the closed-form variance of the weighting and augmented weighting estimators from the sample variance of the estimated influence function. Furthermore, we showed that the influence function of the weighting estimators can be regarded as the influence function of the unadjusted estimator minus its projection to the nuisance tangent space, and proved that the asymptotic variance of the propensity score weighting estimators is no larger than that of the unadjusted estimator. This has important implications, as it provides a compelling rationale for covariate adjustment when targeting win estimands in RCTs, reinforcing the principle that baseline covariates should almost always be leveraged when available rather than ignored. More broadly, this finding strengthens the position reflected in the FDA guidance for COVID-19 trials (U.S. Food and Drug Administration, 2023b), which stated: “To improve the precision of treatment effect estimation and inference, sponsors should consider adjusting for prespecified prognostic baseline covariates …and should propose methods of covariate adjustment.” Our work thus illustrates the operational steps involved in estimating win estimands with the unadjusted and covariate-adjusted win-statistic estimators in the winPSW R package.

Through extensive simulations, we found that, compared to the unadjusted estimator, propensity score weighting estimators are clearly more statistically efficient by leveraging prognostic baseline covariates for estimating both the WR and WD. Furthermore, the augmented weighting estimators showed a further improvement in efficiency, except under unbalanced randomization and with correctly-specified outcome regression based on the interaction model (5.2); in that case, the performance of the propensity score weighting and augmented weighting estimators performed rather similarly, even when the sample size increases. When the sample size was relatively small, e.g., $n\leq 120$, the OW estimator is marginally more efficient than the IPW estimator. Under an unbalanced design with smaller sample sizes, AOW is also more efficient than AIPW and empirical type I error rate of AOW is closer to the nominal 5% level. But in most cases with sufficient sample sizes, our observations are that OW and IPW, AOW and AIPW have similar performance when addressing the WR and WD, which is consistent with prior findings for comparing IPW and OW in randomized trials with a binary outcome (Zeng et al., 2021). Finally, as a cautionary note, we observed that when estimating the WR, the IPW estimator could be less efficient than the unadjusted estimator when the sample size was extremely small (e.g., $n=60$) and the degree of treatment effect heterogeneity was moderate. Based on this observation, we advocate using the OW estimator when the sample size is limited, to safeguard potential finite-sample efficiency loss.

Two estimands, the WR and WD, served as our running examples for demonstrating efficiency gains from propensity score weighting and augmented weighting. We caution, however, that the WR alone may be a misleading summary of an ordinal treatment effect in the presence of ties, and this limitation has been the subject of active debate (Ajufo et al., 2023; Butler et al., 2024; Davison et al., 2025). This is because WR is the ratio of winning to losing comparisons, and is invariant to the absolute number of decisive comparisons and therefore does not always accurately reflect the magnitude of evidence. For example, a trial with 3 wins, 1 loss, and 4996 ties and a trial with 3000 wins, 1000 losses, and 1000 ties both yield a WR of 3, even though the latter conveys a far stronger treatment benefit. The WD does not share this insensitivity, as it normalizes by all pairwise comparisons rather than only decisive ones, and would distinguish the two trials above (0.0004 versus 0.40). For this reason, we recommend reporting the WR and WD together when analyzing ordinal outcomes, rather than relying on the WR alone, so that both the relative and absolute views of the treatment effect are conveyed. Alternatively, a useful measure that incorporates ties is the win odds (WO, or success odds) (Dong et al., 2020b; Brunner et al., 2021), defined as

[TABLE]

After incorporating tied comparisons into the analysis, the WO not only provides a more accurate reflection of similarity between groups and ensures greater numerical stability, but also improves interpretation since it can be interpreted as a true odds (Verbeeck et al., 2023). In principle, our covariate-adjusted estimators can be generalized to accommodate WO. To proceed, the win probability defined in (2.1) should first be modified to be the probabilistic index, that is, $\tau_{1}=\mathbb{P}\left\{Y_{i}(1)>Y_{j}(0)\right\}+\mathbb{P}\left\{Y_{i}(1)=Y_{j}(0)\right\}/2$. Second, when augmented weighting estimators are used, the outcome model $\mu(\cdot,\cdot;\alpha)$ could be fitted by the more general, probabilistic index model (PIM) proposed by Thas et al. (2012). Then, by using the same procedure for weighting and augmented weighting estimators described in Section 3 and Section 4, we can obtain corresponding WO estimators. We expect similar patterns for the covariate-adjusted WO estimator but defer a full development to future research. Finally, Harrell (2021) proposed a Markov proportional odds ordinal logistic model that analyzes longitudinal ordinal outcomes conditional on baseline covariates (including treatment) and the previous state, allowing the treatment effect to be estimated as a function of time. There are, however, two key distinctions from our work. Above all, their approach requires the ordinal outcome to be measured repeatedly over time for each individual, whereas our estimators target a single cross-sectional ordinal outcome. In addition, the native parameter of the transition model is a within-patient odds ratio conditional on the previous state, which, even after marginalization, yields time-indexed state-occupancy probabilities rather than the marginal, between-patient pairwise-comparison estimands considered in our framework. Although their model can provide more efficient derived estimates than simple patient-level summary statistics under a correctly specified transition model, the differing estimands and data structures mean that a direct efficiency comparison is not straightforward, and we leave a formal investigation to future research.

We conclude by summarizing two promising future research directions. First, although our study focuses on the weighting and augmented weighting estimators in the context of randomized trials, the proposed estimators and software are equally applicable to observational studies with ignorable treatment assignment, though the target estimand based on pairwise comparisons will differ depending on the choice of the tilting function $h(X_{i},X_{j})$; see Li et al. (2018) for a detailed elaboration in the context of the weighted average treatment effect estimand (not constructed based on pairwise comparisons). For example, although IPW targets the entire combined population from treatment and control, OW would target the subpopulation with sufficient covariate overlap and clinical equipoise (Li et al., 2018). By shifting the goalpost, OW has theoretically and empirically proven to be the most efficient choice among the class of balancing weights (Li et al., 2018; Li and Li, 2019; Cao et al., 2025), and is the recommended approach to address lack of positivity in the treatment assignment (Li et al., 2019). However, the optimality of OW has so far been established only for estimating the average treatment effect; additional work is needed to establish its optimality for estimating pairwise comparison estimands in the form of (3.1). Second, although we have focused on a univariate ordinal outcome, the win-statistics framework is most often applied to prioritized composite endpoints (Buyse et al., 2025), and our estimators can be extended to this setting. That is, as long as the win and loss functions are well defined, the machinery of Sections 3 and 4 carries over directly, with the augmented estimators additionally requiring a plausible working outcome model (e.g., through the PIM). A key complication arises when the prioritized components involve right-censored time-to-event outcomes, where the win and loss estimands are defined in the absence of censoring. In such cases, the inverse probability of censoring weighting has been suggested in the literature (Dong et al., 2020a; Cao et al., 2026; Fang and Li, 2026), and a full development of covariate-adjusted win statistics for censored composite endpoints is left to future research.

Acknowledgement

F. Li, G. Tong and PJ Heagerty are supported by the United States National Institutes of Health (NIH), National Heart, Lung, and Blood Institute (NHLBI, grant number 1R01HL178513). K. Davis-Plourde is supported by CTSA Grant Number UL1 TR001863 from the National Center for Advancing Translational Science (NCATS), a component of the National Institutes of Health (NIH). All statements in this report, including its findings and conclusions, are solely those of the authors and do not necessarily represent the views of the NIH.

Data Availability Statement

The de-identified data from the ORCHID study are publicly available from the NHLBI BioLINCC website and can be requested at https://biolincc.nhlbi.nih.gov/studies/petal_orchid/. F. Li obtained access to the publicly available ORCHID study data according to a Data Use Agreement from BioLINCC and performed the data analysis in Section 6.

Supporting information

R code used in this article for simulation studies can be found at https://github.com/Zhiqiangcao/winPSW/tree/main/simulation_code. Additional supporting information including Web Appendices A–F can be found at Github platform https://github.com/Zhiqiangcao/winPSW.