Computation of The Mann-Whitney Effect under Parametric Survival Copula Models

1. Introduction

When comparing the survival times of two independent groups, the Mann-Whitney parameter plays an important role in the two-sample problem ([1]). The Mann-Whitney parameter, say p, is defined as the probability that a random subject from one group (with survival time $T_1$ in group 1) survives longer than an independent random subject from the other group (with survival time $T_2$ in group 2), plus one-half the probability that the two subjects survive at the same time:

$$p=\mathbb{P}(T_1>T_2)+\frac{1}{2}\mathbb{P}(T_1=T_2).$$

(1)

The Mann-Whitney effect relates to important statistical ideas, such as, the Mann-Whitney test ([2]), hazard ratios, and win ratio ([3]). The Mann-Whitney test examines the null hypothesis $H_0:p=1/2$ v.s. $H_1:p\ne 1/2$. The hazard ratio is the main effect measure of a Cox proportional hazards model, which is a typical statistical model in survival analysis. The win ratio w is given by the odds of p; that is, $w=p/(1-p)$. That is, $p>1/2$, or equivalently $w>1$, implies a protective survival effect for group 1.

The problem of estimating the parameter p plays an important part in survival analysis. The basic idea was first studied by [4]. They illustrated an attractive relationship between the Mann-Whitney statistic and the stress-strength model. [1] first proposed a nonparametric estimator for p under independent censoring. Since then, this topic has been investigated by several researchers. In the following, we refer to some recent studies in the field of survival analysis. [5] modified Efron’s estimator for p under small sample sizes. [6] proposed a copula-graphic estimator for p and suggested the Mann-Whitney test to compare two survival distributions in the presence of dependent censoring. [7] introduced the Bayesian estimation of p for the log-Lindley distribution. [8] proposed estimating the Mann-Whitney effects in factorial clustered data. [9] developed methodologies for constructing fixed-accuracy confidence intervals of p when $T_1$ and $T_2$ follow geometric and the exponential distributions, respectively. [10] studied a estimation procedure of the stress-strength model for the two independent unit-half-normal distributions with different shape parameters. [11] investigated the effect of dependence of the valiables on p in the stress-strength model with the exponential margins. [12] proposed a group sequential method for estimating the Mann-Whitney parameter. [13] studied the estimator of p in point, interval, and Bayesian estimations when the stress valiables follow geometric and Lindley distribution. All the methods assumed that $T_1$ and $T_2$ are independent.

When $T_1$ and $T_2$ are independent of each other and continuous, one can estimate p with the marginal distributions based on the following integral:

$$p=-{\int }_0^{\infty }{S}_1\left(t\right)d{S}_2\left(t\right).$$

That is, one can estimate p by estimating two marginal survival functions $S_1$ and $S_2$. However, this is not the case when $T_1$ and $T_2$ are dependent; the phenomenon is sometimes called “Hand’s paradox” ([14]). This showed that the paradox arises when $T_1$ and $T_2$ are regarded as potential outcomes in the framework of causal inference. Therefore, p in the integral cannot be interpreted as the true treatment effect. Besides, dependence of outcomes from observation to observation is well-known in factorial designs and cross-over designs ([15]).

Since p is not identifiable solely from independently sampled data, [16] suggested a bound for p under all possible dependence structures for $T_1$ and $T_2$. Alternatively, [17] reformulated p such that it can be identified from randomized treatment assignments. However, to estimate the true p, we must model the bivariate survival function of $T_1$ and $T_2$. Copula is often used to model joint distributions of dependent survival times ([18,19]).

In this article, we propose a model for the bivariate survival function by using parametric copulas and parametric marginal distributions. We then derive a new formula for computing p by a one-dimensional integral. We also propose a new formula for p under the restricted follow-up. To make the proposed computation method for p to be easily performed by users, we develop a Shiny-based web app. Furthermore, we validate the accuracy of the proposed computation method and Shiny web app by simulations. We finally illustrate the proposed method by two real datasets.

The rest of the paper is organized as follows. In Section 2, we review copula-based models and introduce several well-known copula families. In this section, we show that one can compute p by Theorem in [20], and we extend the theorem to compute p when the follow-up time is restricted up to time $\tau$. In Section 3, we introduce a Shiny web app in the R that can compute p via simple commands. In Section 4, we describe a simulation study to show the correctness of the proposed calculator for p. In Section 5, we illustrate a meaningful application of our proposed method using survival data.

2. Proposed Method

In this section, we first introduce copula-based bivariate survival models for $T_1$ and $T_2$. We then propose our method for computing the Mann-Whitney effect p in Equation (1) under the copula models.

2.1. Copula Models for Dependent Survival Time

According to [21], any bivariate distribution function for $(T_1,T_2)$ can be formulated by using a copula. A bivariate copula is a bivariate distribution function for two uniform variables on $[0,1]$ ([22,23]). Let $T_1$ and $T_2$ be continuous survival times with marginal survival functions $S_1$ and $S_2$, respectively. Initially, we model the bivariate survival function of $T_1$ and $T_2$ using a copula C:

$$\mathbb{P}(T_1>t_1,T_2>t_2)=C(S_1\left(t_1\right),S_2\left(t_2\right)).$$

(2)

This representation is useful since two marginal survival functions $S_1$ and $S_2$ are separated from the dependence structure C.

We will consider the following well-known families of bivariate copulas.

• The independence copula:

  $$C(u,v)=uv.$$
• The Clayton copula ([24]):

  $$C_{\theta }(u,v)=max\left({(u^{-\theta }+v^{-\theta }-1)}^{-1/\theta },0\right),\theta \in [-1,\infty )\setminus \left\{0\right\}.$$
• The Gumbel copula ([25]):

  $$C_{\theta }(u,v)=exp\left\{-{[{(-logu)}^{\theta +1}+{(-logv)}^{\theta +1}]}^{\frac{1}{\theta +1}}\right\},\theta \in [0,\infty ).$$
• The Frank copula ([26]):

  $$C_{\theta }(u,v)=-\frac{1}{\theta }log\left(1+\frac{(e^{-\theta u}-1)(e^{-\theta v}-1)}{e^{-\theta }-1}\right),\theta \in (-\infty ,\infty )\setminus \left\{0\right\}.$$
• The Farlie-Gumbel-Morgenstern (FGM) copula ([27]):

  $$C(u,v)=uv+\theta uv(1-u)(1-v),\theta \in [-1,1].$$
• The Gumbel-Barnett (GB) copula ([22,28,29]):

  $$C(u,v)=uvexp(-\theta logulogv),\theta \in [0,1].$$

In Figure 1, we present scatterplots generated from various copulas with parameter $\theta$. The Clayton copula (Figure 1 (b), (c)) shows lower tail dependence, the Gumbel copula (Figure 1 (d)) shows upper tail dependence, Frank copula (Figure 1 (e), (f)) shows symmetric dependence around the median, and the FGM copula (Figure 1 (g), (h)) is similar to the Frank copula. Unlike other copulas, the GB copula exhibits negative dependence only (Figure 1 (i)).

These copulas have been applied to survival data and other data analyses. The Clayton copula was applied to survival data with dependent censoring. For instance, [30] modeled dependence between survival and dependent censoring times in the survival data of tuberculosis cure. The Clayton, Gumbel, and Frank copula were also applied to dependently censored data in clinical trials or observational studies ([31,32,33]). The Gumbel, Frank, and FGM copulas were often used in competing risks models on survival data analysis ([34,35,36]). Copulas were also applied to multivariate meta-analysis; [37] and [38] proposed a bivariate Clayton, FGM and Gumbel models for bivariate meta-analysis. The Gumbel-Barnett copula has the simple form and is suitable for modeling negative dependence ([22,28]). Therefore, it is important to consider a variety of copulas for dependent survival times.

To see the strength of dependence in a copula, $\theta$ can be transformed to Kendall’s tau, $\tau$. Kendall’s tau is a well-known measure to assess the dependence between two variables. Under the copula model (2), Kendall’s tau for $T_1$ and $T_2$ is expressed as

$$\mathrm{Kendall}\text{’}\mathrm{s}\tau =4{\int }_0^1{\int }_0^1{C}_{\theta }(u,v)C_{\theta }(du,dv)-1.$$

Kendall’s $\tau$ does not depend on the marginals, and is solely determined by the copula. Therefore, it is advantageous over the Pearson correlation for $T_1$ and $T_2$. Kendall’s $\tau$ of each copula is expressed as follow.

• The independence copula:

  $$\mathrm{Kendall}\text{’}\mathrm{s}\tau =0.$$
• The Clayton copula:

  $$\mathrm{Kendall}\text{’}\mathrm{s}\tau =\frac{\theta }{\theta +2},\theta \in [-1,\infty )\setminus \left\{0\right\}.$$
• The Gumbel copula:

  $$\mathrm{Kendall}\text{’}\mathrm{s}\tau =\frac{\theta }{\theta +1},\theta \in [0,\infty ).$$
• The Frank copula:

  $$\mathrm{Kendall}\text{’}\mathrm{s}\tau =1-\frac{4}{\theta }+\frac{4D_{\theta }}{\theta },\mathrm{where}{D}_{\theta }=\frac{1}{\theta }{\int }_0^{\theta }\frac{x}{exp\left(x\right)-1}dx,\theta \in (-\infty ,\infty )\setminus \left\{0\right\}.$$
• The FGM copula:

  $$\mathrm{Kendall}\text{’}\mathrm{s}\tau =\frac{2}{9}\theta ,\theta \in [-1,1].$$
• The GB copula:

  $$\mathrm{Kendall}\text{’}\mathrm{s}\tau =1-\frac{4}{\theta }{\int }_0^{1}t(1-\theta logt)log(1-\theta logt)dt,\theta \in [0,1].$$

2.2. Proposed Method for Computing p

In this section, we propose a new formula for computing p under the survival copula model (2). Let $U=S_1\left(T_1\right)$ and $V=S_2\left(T_2\right)$. For computing p, we will use the conditional distribution function for U given $V=v$, which is the partial derivative of C with respect to v:

$$C^{[0,1]}(u,v)=\mathbb{P}(U\le u\mid V=v)=\frac{\partial C(u,v)}{\partial v}.$$

Then, by slightly modifying Theorem of [20], p can be expressed as the univariate integral on $[0,1]$:

$$\begin{array}{cc}\hfill p& =\mathbb{P}(T_1>T_2)+\frac{1}{2}\mathbb{P}(T_1=T_2)\hfill \\ \hfill & =\mathbb{P}(S_1\left(T_1\right)<S_1\left(T_2\right))=\mathbb{P}(U<S_1\left(S_2^{-1}\left(V\right)\right))=\mathbb{E}\left[\mathbb{P}(U<S_1\left(S_2^{-1}\left(V\right)\right)\mid V)\right]\hfill \\ \hfill & ={\int }_0^1{C}^{[0,1]}(S_1\left(S_2^{-1}\left(v\right)\right),v)dv.\hfill \end{array}$$

(3)

We note that Theorem of [20] is not directly applicable to the survival copula model (2) since that theorem is designed for the copula model $\mathbb{P}(T_1\le t_1,T_2\le t_2)=C(1-S_1\left(t_1\right),1-S_2\left(t_2\right))$.

Equation (3) with aforementioned copulas is computed by the following formulas:

• The Clayton copula

  $$p={\int }_0^1{v}^{-\theta -1}{\left({\left\{S_1\left(S_2^{-1}\left(v\right)\right)\right\}}^{-\theta }+v^{-\theta }-1\right)}^{-\frac{1}{\theta }-1}dv.$$
• The Gumbel copula

  $$\begin{array}{cc}\hfill p=& {\int }_0^1\left\{exp\left\{-{[{(-log{S}_1\left(S_2^{-1}\left(v\right)\right))}^{\theta +1}+{(-logv)}^{\theta +1}]}^{\frac{1}{\theta +1}}\right\}\right.\hfill \\ \hfill & \times \left.{\left[{\left(-log\left(S_1\left(S_2^{-1}\left(v\right)\right)\right)\right)}^{\theta +1}+{\left(-log\left(v\right)\right)}^{\theta +1}\right]}^{-\frac{\theta }{\theta +1}}\frac{{(-log\left(v\right))}^{\theta }}{v}\right\}dv.\hfill \end{array}$$
• The Frank copula

  $$p={\int }_0^1\frac{e^{-\theta v}\left(e^{-\theta S_1\left(S_2^{-1}\left(v\right)\right)}-1\right)}{\left(e^{-\theta }-1\right)+\left(e^{-\theta S_1\left(S_2^{-1}\left(v\right)\right)}-1\right)\left(e^{-\theta v}-1\right)}dv.$$
• The FGM copula

  $$p={\int }_0^1\left\{S_1\left(S_2^{-1}\left(v\right)\right)+\theta S_1\left(S_2^{-1}\left(v\right)\right)\left(1-S_1\left(S_2^{-1}\left(v\right)\right)\right)(1-2v)\right\}dv.$$
• The GB copula

  $$p={\int }_0^1\left\{S_1\left(S_2^{-1}\left(v\right)\right)\left(1-\theta log{S}_1\left(S_2^{-1}\left(v\right)\right)\right)v^{-\theta log{S}_1\left(S_2^{-1}\left(v\right)\right)}\right\}dv.$$

In order to compute p by the above formulas, we need to specify $S_1$, $S_2$, and $\theta$. One can specify $S_1$ and $S_2$ by continuous parametric models that will be discussed in Section 2.4. One can try different values for $\theta$ in a sensitivity analysis. Note that the above calculations are not applicable for discrete parametric models for $S_1$ and $S_2$.

2.3. Computing p with Follow-Up Time

For survival data, the follow-up period is often limitted. When a subject survives longer than the follow-up period, one may treat the survival time of the subject as equal to the follow-up period ([5,39,40]). In this section, we assume that every subject has a common follow-up time $\tau$. This means that we define the Mann-Whitney effect for $min(T_1,\tau )$ and $min(T_2,\tau )$. We now obtain p with follow-up time $\tau$ from the following theorem, a straightforward expansion of Theorem in [20].

Theorem 1.

The Mann-Whitney effect p with a follow-up time τ is written as the univariate integral:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill p_{\tau }& =\mathbb{P}(min(T_1,\tau )>min(T_2,\tau ))+\frac{1}{2}\mathbb{P}(min(T_1,\tau )=min(T_2,\tau ))\hfill \\ \hfill & =\mathbb{P}(T_1>T_2,T_2<\tau )+\frac{1}{2}\mathbb{P}(T_1>\tau ,T_2>\tau )\hfill \\ \hfill & =\mathbb{P}(U<S_1\left(S_2^{-1}\left(V\right)\right),V>S_2\left(\tau \right))+\frac{1}{2}C(S_1\left(\tau \right),S_2\left(\tau \right))\hfill \\ \hfill & ={\int }_{S_2\left(\tau \right)}^1{C}^{[0,1]}(S_1\left(S_2^{-1}\left(v\right)\right),v)dv+\frac{1}{2}C(S_1\left(\tau \right),S_2\left(\tau \right)).\hfill \end{array}\end{array}$$

(4)

Equation 4 with different copulas is computed by the following formulas:

• The Clayton copula

  $$p_{\tau }={\int }_{S_2\left(\tau \right)}^1\left\{v^{-\theta -1}{\left({\left\{S_1\left(S_2^{-1}\left(v\right)\right)\right\}}^{-\theta }+v^{-\theta }-1\right)}^{-\frac{1}{\theta }-1}\right\}dv+\frac{1}{2}{(S_1{\left(\tau \right)}^{-\theta }+S_2{\left(\tau \right)}^{-\theta }-1)}^{-1/\theta }.$$
• The Gumbel copula

  $$\begin{array}{cc}\hfill p_{\tau }=& {\int }_{S_2\left(\tau \right)}^1\left\{exp\left\{-{[{(-log{S}_1\left(S_2^{-1}\left(v\right)\right))}^{\theta +1}+{(-logv)}^{\theta +1}]}^{\frac{1}{\theta +1}}\right\}\right.\hfill \\ \hfill & \left.\times {\left[{\left(-log\left(S_1\left(S_2^{-1}\left(v\right)\right)\right)\right)}^{\theta +1}+{\left(-log\left(v\right)\right)}^{\theta +1}\right]}^{-\frac{\theta }{\theta +1}}\frac{{(-log\left(v\right))}^{\theta }}{v}\right\}dv\hfill \\ \hfill & +\frac{1}{2}exp\left\{-{[{(-log{S}_1\left(\tau \right))}^{\theta +1}+{(-log{S}_2\left(\tau \right))}^{\theta +1}]}^{\frac{1}{\theta +1}}\right\}.\hfill \end{array}$$
• The Frank copula

  $$\begin{array}{cc}\hfill p_{\tau }=& {\int }_{S_2\left(\tau \right)}^1\frac{e^{-\theta v}\left(e^{-\theta S_1\left(S_2^{-1}\left(v\right)\right)}-1\right)}{\left(e^{-\theta }-1\right)+\left(e^{-\theta S_1\left(S_2^{-1}\left(v\right)\right)}-1\right)\left(e^{-\theta v}-1\right)}dv\hfill \\ \hfill & +\frac{1}{2}\left(-\frac{1}{\theta }log\left(1+\frac{(e^{-\theta S_1\left(\tau \right)}-1)(e^{-\theta S_2\left(\tau \right)}-1)}{e^{-\theta }-1}\right)\right).\hfill \end{array}$$
• The FGM copula

  $$\begin{array}{cc}\hfill p_{\tau }=& {\int }_{S_2\left(\tau \right)}^1\left\{S_1\left(S_2^{-1}\left(v\right)\right)+\theta S_1\left(S_2^{-1}\left(v\right)\right)\left(1-S_1\left(S_2^{-1}\left(v\right)\right)\right)(1-2v)\right\}dv\hfill \\ \hfill & +\frac{1}{2}\left(S_1\left(\tau \right)S_2\left(\tau \right)+\theta S_1\left(\tau \right)S_2\left(\tau \right)(1-S_1\left(\tau \right))(1-S_2\left(\tau \right))\right).\hfill \end{array}$$
• The GB copula

  $$\begin{array}{cc}\hfill p_{\tau }=& {\int }_{S_2\left(\tau \right)}^1\left\{S_1\left(S_2^{-1}\left(v\right)\right)\left(1-\theta log{S}_1\left(S_2^{-1}\left(v\right)\right)\right)v^{-\theta log{S}_1\left(S_2^{-1}\left(v\right)\right)}\right\}dv\hfill \\ \hfill & +\frac{1}{2}\left(S_1\left(\tau \right)S_2\left(\tau \right)exp(-\theta log{S}_1\left(\tau \right)log{S}_2\left(\tau \right))\right).\hfill \end{array}$$

2.4. Marginal Survival Distributions

To compute p and $p_{\tau }$, we considered following three parametric distributions as the marginal survival distributions of group $j\in \{1,2\}$:

• The exponential distribution:

  $$\begin{array}{cc}\hfill & S_j\left(t\right)=exp(-{\lambda }_{j}t),{\lambda }_j>0,\hfill \\ \hfill & S_i\left(S_j^{-1}\left(v\right)\right)=exp\left(\frac{{\lambda }_i}{{\lambda }_j}logv\right).\hfill \end{array}$$

  where ${\lambda }_j$ is a rate parameter.
• The Weibull distribution:

  $$\begin{array}{cc}\hfill & S_j\left(t\right)=exp(-{\lambda }_j{t}^{k_j}),{\lambda }_j>0,k_j>0,\hfill \\ \hfill & S_i\left(S_j^{-1}\left(v\right)\right)=exp\left(-{\lambda }_i{\left(-\frac{logv}{{\lambda }_j}\right)}^{\frac{k_i}{k_j}}\right).\hfill \end{array}$$

  where ${\lambda }_j$ is a scale parameter and $k_j$ is a shape parameter.
• The gamma distribution:

  $$S_j\left(t\right)=1-\frac{\gamma (k_j,{\lambda }_{j}t)}{\Gamma \left(k_j\right)},{\lambda }_j>0,k_j>0.$$

  where ${\lambda }_j$ is a scale parameter and $k_j$ is a shape parameter, and $\Gamma \left(k\right)$ is the gamma function, and $\gamma (k,\lambda t)={\int }_0^{\lambda t}{x}^{k-1}{e}^{-x}dx$ is the lower incomplete gamma function. The gamma distribution has no simple closed-form expression for the inverse survival function. Therefore, one can use approximations for the inverse survival function. In Section 3, we use the R function “qgamma” to calculate the quantile funciton of the gamma distribution.

In Figure 2, we present survival curves of three parametric distributions, the exponential, Weibull, and gamma distributions with different parameters. These plots show that these distributions can represent almost any continuous survival curve that will be encountered in practice.

Example 1.

Let the marginals $S_1\left(t\right),S_2\left(t\right)$ be the exponential distributions with parameter ${\lambda }_1=1,{\lambda }_2=2$ and $T_1$ and $T_2$ are independent. Then by Theorem 1 with $\tau =\infty$, p is given by

$$\begin{array}{cc}\hfill p& ={\int }_0^1{C}^{[0,1]}(S_1\left(S_2^{-1}\left(v\right)\right),v)dv\hfill \\ \hfill & ={\int }_0^{1}exp\left(\frac{{\lambda }_1}{{\lambda }_2}logv\right)dv=\frac{{\lambda }_2}{{\lambda }_1+{\lambda }_2}=\frac{2}{3}.\hfill \end{array}$$

□

Example 2.

Let the marginals $S_1\left(t\right),S_2\left(t\right)$ be the exponential distributions with parameter ${\lambda }_1=1,{\lambda }_2=2$ and $C(u,v)$ be the Clayton copula with parameter $\theta =3$. Then Kendall’s τ is given by

$$\mathrm{Kendall}\text{’}\mathrm{s}\tau =\frac{\theta }{\theta +2}=0.6$$

and by Theorem 1 with $\tau =\infty$, p is given by

$$\begin{array}{cc}\hfill p& ={\int }_0^1{C}^{[0,1]}(S_1\left(S_2^{-1}\left(v\right)\right),v)dv\hfill \\ \hfill & ={\int }_0^1{v}^{-\theta -1}{\left(v^{-\theta \frac{{\lambda }_1}{{\lambda }_2}}+v^{-\theta }-1\right)}^{-\frac{1}{\theta }-1}dv\hfill \\ \hfill & ={\int }_0^1{v}^{-4}{\left(v^{-3\frac{{\lambda }_1}{{\lambda }_2}}+v^{-3}-1\right)}^{-\frac{4}{3}}dv=0.84.\hfill \end{array}$$

When $\tau <\infty$, $p_{\tau }$ is given by

$$\begin{array}{cc}\hfill p_{\tau }& ={\int }_0^1{C}^{[0,1]}(S_1\left(S_2^{-1}\left(v\right)\right),v)dv+\frac{1}{2}\mathbb{P}(min(T_1,\tau )=min(T_2,\tau ))\hfill \\ \hfill & ={\int }_{e^{2\tau }}^1{v}^{-4}{\left(v^{-3\frac{{\lambda }_1}{{\lambda }_2}}+v^{-3}-1\right)}^{-\frac{4}{3}}dv+{\left(e^{-3\tau }+e^{-6\tau }-1\right)}^{-\frac{1}{3}}.\hfill \end{array}$$

When $\tau =0.5$, $p_{\tau }=0.68$. Here, we computed the last two equations by a numerical integration by the R function "integrate". □

Parameters for the marginal distributions can be estimated by maximum likelihood estimators (MLEs) when the survival data are available in two groups (Section 5).

3. Software and Web App

We developed a Shiny-based web app to implement the proposed method of computing p. The Shiny app is available at (https://nkosuke.shinyapps.io/shiny_survival/) and can be used by any environment, including smartphones. Using this web app, users can choose a marginal survival distribution, copula, and the relevant parameters to compute p. Our web app is easy to use without knowledge of the R.

3.1. Input

We considered three marginal survival distributions; the exponential, Weibull, and gamma distributions. We considered the Clayton, Gumbel, Frank, FGM, and GB copulas. One can choose marginal distributions and copulas, set their parameters, and choose the langage displayed on the screen in the input panels on the left-hand of our app (Figure 3). Furthermore, one can set a follow-up time $\tau$ with a slider bar. In Figure 3, we set the marginal survival distribution = “Exponential Distribution”, ${\lambda }_1=0.5$, ${\lambda }_2=0.25$, $\tau =4.5$, copula = “Clayton”, $\theta =1.5$ and the langage displayed on the screen “English”.

3.2. Output

This app displays several formulas, survival curves of each group, p, $p_{\tau }$, and Kendall’s $\tau$. These formulas include marginal and bivariate survival functions and formulas of p and $p_{\tau }$ based on the input values. The theoretical values of p and $p_{\tau }$ are displayed together with survival curves.

Figure 3 displays the web app. In this settings, we have the output $p=0.22,p_{\tau }=0.27$, Kendall’s $\tau =0.429$.

4. Simulation Studies

To show the correctness of the proposed calculator for p, we conduced a simulation study. For the simulation study, we set the marginal survival functions to be the exponential distributions with $({\lambda }_1,{\lambda }_2)=(1,2)$, the Weibull distributions with $({\lambda }_1,k_1,{\lambda }_2,k_2)=(1,0.5,2,1)$, and the gamma distributions with $({\lambda }_1,k_1,{\lambda }_2,k_2)=(1,1.5,2,2)$. Furthermore, we set the copula parameters, the Clayton copula with $\theta =1,5\mathrm{or}10$, the Gumbel copula with $\theta =0\mathrm{or}4$, the Frank copula with $\theta =-5,1,\mathrm{or}5$, the FGM copula with $\theta =-1,0,\mathrm{or}1$, the GB copula with $\theta =0.5,\mathrm{or}1$, and follow-up time $\tau =0.5,2,\mathrm{or}5$. We generated 100,000 pairs $(T_{i1},T_{i2}),i=1,\dots ,M,\mathrm{where}M=100,000$, from the bivariate survival function based on the aforementioned settings, and calculated the Monte Carlo simulation values defined as,

$$\begin{array}{c}\hfill p_{\tau ,\mathrm{sim}}=\frac{1}{M}\sum _{i=1}^M\left[\mathbb{1}(min(T_{i1},\tau )>min(T_{i2},\tau ))+\frac{1}{2}\mathbb{1}(min(T_{i1},\tau )=min(T_{i2},\tau ))\right].\end{array}$$

Table 1 shows that the simulation values are nearly equal to the theoretical values computed by the formula of Theorem 1 for every setting. In conclusion, our simulations show that Theorem 1 is correct, and the Shiny web app based on Theorem 1 is reliable.

5. Numerical Examples

In this section, we apply our proposed methods to a tongue cancer dataset and a prostate cancer dataset. Before analyzing the real datasets, we introduce basic notations and ideas for estimating p by using censored data. Let $T_{ij},j=1,2,i=1,\dots ,n_j$ be survival times, $C_{ij},j=1,2$ be censored times, $t_{ij}=min(T_{ij},C_{ij}),j=1,2$ be observed time, and ${\Delta }_{ij}=\mathbb{1}(T_{ij}\le C_{ij}),j=1,2$ be event indicator. What we observe is $t_{ij}$ and ${\Delta }_{ij}$. That is, ${\Delta }_{ij}$ is 0 or 1 according to whether $t_{ij}$ is a censored time or a survival time. As the exponential distribution is shown to fit well for $T_{i1}$ and $T_{i2}$, we obtained MLE of the exponential hazard rate ${\widehat{\lambda }}_j,j=1,2$ by

$${\widehat{\lambda }}_j=\frac{{\sum }_{i=1}^{n_j}{\Delta }_{ij}}{{\sum }_{i=1}^{n_j}{t}_{ij}},j=1,2.$$

Then, by applying the values of the MLE to the proposed Shiny web app, we obtained the estimators $\widehat{p}$ and ${\widehat{p}}_{\tau }$, where $\tau$ was chosen appropriately (Section 5.1 and Section 5.2). On the other hand, under the independence assumption of $T_{i1}$ and $T_{i2}$, the naïve estimator of p is

$$\begin{array}{cc}\hfill {\widehat{p}}^{\mathrm{KM}}& =-{\int }_0^{\infty }{\widehat{S}}_1^{\pm }\left(t\right)d{\widehat{S}}_2\left(t\right),\hfill \\ \hfill {\widehat{p}}_{\tau }^{\mathrm{KM}}& =-{\int }_0^{{\widehat{S}}_2^{-1}\left(\tau \right)}{\widehat{S}}_1^{\pm }\left(t\right)d{\widehat{S}}_2\left(t\right)+\frac{1}{2}{\widehat{S}}_1^{\pm }\left(\tau \right){\widehat{S}}_2^{\pm }\left(\tau \right),\hfill \end{array}$$

where ${\widehat{S}}_i^{\pm }\left(t\right)=[{\widehat{S}}_i(t+)+{\widehat{S}}_i(t-)]/2$ and ${\widehat{S}}_i\left(t\right)$ is a Kaplan-Meier (KM) estimator. However, this estimate is subject to the independence of two groups. Therefore, the proposed estimator is useful to examine the sensitivity under a variety of dependence structures via copulas.

5.1. Tongue Cancer Data

The tongue dataset is available in the R package KMsurv. It has 80 observations and contains: type (Tumor DNA profile: 1 = aneuploid tumor, 2 = diploid tumor), time (Time to death or on-study time (weeks)), and death (Event indicator: 0 = alive, 1 = dead). It contains $n_1=52$ observations in aneuploid cancer group ($j=1$), and $n_2=28$ observations in diploid cancer group ($j=2$). We considered the follow-up time $\tau ={min}_i\left\{{max}_j{t}_{ij}{\Delta }_{ij}\right\}$ and obtained $\tau =167$. The tongue cancer data resulted in ${\widehat{\lambda }}_1=0.00736,{\widehat{\lambda }}_2=0.0130$ and ${\widehat{p}}_{\tau }^{\mathrm{KM}}=0.615$. In Figure 4, the KM estimators of each groups and the estimated exponential survival curves are plotted. We conducted sensitivity analyses using copula-based approach. We calculated $\widehat{p}$ by Theorem 1 under weak, strong positive, and negative independences. We calculated ${\widehat{p}}_{\tau }$ via the web app (Section 3). Figure 5 shows the output under the independent, Clayton, Gumbel, Frank, FGM, and GB copula with parameter $\theta \in \{1,5\}\left(\mathrm{the}\mathrm{Clayton}\right),\theta =4\left(\mathrm{the}\mathrm{Gumbel}\right),\theta \in \{-5,5\}\left(\mathrm{the}\mathrm{Frank}\right),\theta \in \{-1,1\}\left(\mathrm{the}\mathrm{FGM}\right),\theta \in \{0.5,1\}\left(\mathrm{the}\mathrm{GB}\right)$. The results under all copulas are summarized in Table 2. We obtained the ${\widehat{p}}_{\tau }$ ranged from $0.600$ to $0.862$ and concluded that a subject in DNA-aneuploid tumor gruop survives longer than in DNA-diploid tumor group. This conclusion did not change under any depencence structures we conducted.

5.2. Prostate Cancer Data

The prostate cancer data is avalible in the R package asaur ([41]). It has 14,294 observations and contains: grade (moderately differentiated and poorly differentiated), survTime (time from diagnosis to death or last date known alive), and status (Event indicator: 0 = censored, 1 = death from prostate cancer). It contains $n_1=\mathrm{10,988}$ observations in moderately differentiated group ($j=1$), and $n_2=\mathrm{3,306}$ observations in poorly differentiated group ($j=2$).

The prostate cancer data resulted in ${\widehat{\lambda }}_1=0.000817,{\widehat{\lambda }}_2=0.00374$, $\tau =108$, and $\widehat{p}=0.679$. In Figure 6, we plot the KM estimators of each group and estimated exponential survival curves. We calculated ${\widehat{p}}_{\tau }^{\mathrm{KM}}$ by Theorem 1 with copulas and several parameters. We calculated ${\widehat{p}}_{\tau }$ under the independent, Clayton, Gumbel, Frank, FGM, and GB copulas with parameter $\theta \in \{1,5\}\left(\mathrm{the}\mathrm{Clayton}\right),\theta =4\left(\mathrm{the}\mathrm{Gumbel}\right),\theta \in \{-5,5\}\left(\mathrm{the}\mathrm{Frank}\right),\theta \in \{-1,1\}\left(\mathrm{the}\mathrm{FGM}\right),\theta \in \{0.5,1\}\left(\mathrm{the}\mathrm{GB}\right)$ via the web app (Figure 7). The results for all scenarios are summarized in Table 3. We obtained the ${\widehat{p}}_{\tau }$ ranged from $0.623$ to $0.665$ and concluded that a subject in moderately differentiated gruop survives longer than in poorly differentiated group. The range of ${\widehat{p}}_{\tau }$ is narrower than one of tongue cancer dataset. The results may be caused by the distinct difference of survival curves between two groups, the short follow-up time for two groups and large sample size.

6. Conclusion

The Mann-Whitney effect has been widely used for survival analysis, which can provide the meaningful measure for treatment effects for survival outcomes. However, the Mann-Whitney effect may not be interpreted as the true treatment effect under dependence of two survival times. In this article, we proposed a parametric copula-based approach for estimating the Mann-Witney effect p under depencence structures for two survival times. We derived the formulas of p under a variety of copulas and marginal survival functions. We also introduced a web-based calculator for p for users. Simulation studies demonstrated the correctness of the proposed calculator for p under a variety of the parametric marginal survival distributions and copulas. The results of data analyses show that the proposed method gives possible changes of p under various denpendence and enables to examine the sensitivity.

In the examples of real datasets, we obtained $p_{\tau }$ under the Clayton, Gumbel, Frank, FGM, GB copulas with varying parameters. The value of $p_{\tau }$ ranged from $0.600$ to $0.862$ in tongue cancer dataset, from $0.623$ to $0.665$ in prostate cancer dataset. We obtained the narrow ranges whose lower bound did not include the null value of $1/2$. The result is consistent with previous studies that Hand’s paradox does not occur under strictly monotonic effect ([14,42]). While more complex dependence structures with various copulas might be considered, the conclusion may not change much.

The main limitation of the present article is that we only discussed the “parametric” approach. However, in practice, researchers may use the “semi-” or “non-parametric” approach. In future work, we will examine the method of computing p without parametric assumptions. Another extension is to include covariates or secondary outcomes in the model, which help obtain narrow bounds for treatment effects ([43]). Another limitation is that only one-parameter copulas are implemented. There are multi-parameter copulas that deserve attention ([22,28,44]).