Factorial survival analysis for treatment effects under dependent censoring

Abstract

Factorial analyses offer a powerful nonparametric means to detect main or interaction effects among multiple treatments. For survival outcomes, for example, from clinical trials, such techniques can be adopted for comparing reasonable quantifications of treatment effects. The key difficulty to solve in survival analysis concerns the proper handling of censoring. So far, all existing factorial analyses for survival data have been developed under the independent censoring assumption, which is too strong for many applications. As a solution, the central aim of this article is to develop new methods for factorial survival analyses under quite general dependent censoring regimes. This will be accomplished by combining existing nonparametric methods for factorial survival analyses with techniques developed for survival copula models. As a result, we will present an appealing F-test that exhibits sound performance in our simulation study. The new methods are illustrated in a real data analysis. We implement the proposed method in an R function surv.factorial(.) in the R package compound.Cox.

Keywords

Copula copula-graphic estimator dependent censoring factorial designs Mann-Whitney test

1 Introduction

Factorial designs allow researchers to assess treatment effects in one-way, two-way, and other layouts of experimental factors.¹ Factorial analyses offer a variety of treatment contrasts involving multiple treatment factors and their levels. The classical analysis of variance (ANOVA) is a test for detecting treatment effects via the F-statistic under normally distributed outcomes. The F-test is an omnibus test^2,3 for the global null hypothesis, where any significance indicates the presence of treatment effects, providing scientific evidence for further explorations of individual factors/treatments.

For survival outcomes, there is a keen interest in adopting factorial analyses for testing treatment effects in cancer clinical trials: see for instance, an ongoing phase III clinical trial for patients with metastatic prostate cancer with a 2 × 2 factorial design.⁴ A 2017 survey identified 30 clinical trials with the 2 × 2 factorial design for phase III cancer treatment,⁵ including a trial that examined the effect of treatments on survival for patients with advanced cervical cancer.⁶

Due to the non-normality of survival data, treatment effects are formulated nonparametrically via the Mann-Whitney effects,⁷ pairwise effects using sample partitions,⁸ the effects on median survival,⁹ or the effects on cumulative hazards.¹⁰ While these metrics provide different survival aspects of treatment effects, they can be identified without the aid of specific distributions. These techniques offer a powerful nonparametric means to detect the main/interaction effects of treatments.

The key difficulty to solve in survival analysis concerns proper handling of censoring owing to the limited follow-up time of clinical trials.¹¹ However, the existing factorial analyses^7–10 were developed under the assumption that the censoring mechanism is independent of survival. Unfortunately, this assumption is imposed for the sake of mathematical convenience, and its validity in real applications is largely questionable.^12–18

Dependent censoring arises in the analysis of survival data, especially in data obtained from biomedical studies.^{14,15,19–29} Typical examples of dependent censoring are the dropout, withdrawal, or removal of patients because of worsening health conditions. For instance, Staplin et al.¹⁵ analyzed the survival of patients registered for an elective liver transplant. Censoring occurs when patients are removed from the waiting list for liver transplants. Transplantation potentially yields dependent censoring because patients closer to death are more likely to receive transplants. Also, patients removed due to deteriorating conditions may have died on the day of removal.

The Kaplan-Meier (KM) estimator for a survival function is systematically biased when censoring is associated with survival time.³⁰ Andersen and Perme¹⁹ suggested reducing the bias of the KM estimator with the aid of observed covariates. Emura and Chen¹⁶ analyzed the bias of Cox regression and variable selection under dependent censoring. Emura and Hsu³¹ studied the bias on two-sample tests based on the Mann-Whitney type treatment effects. Other statistical tools for dependent censoring are reviewed by the books of Emura and Chen¹⁷ and Collett.²⁰ However, dependent censoring has not been studied in the context of factorial designs.

In this article, we propose a new method for factorial survival analyses, which can take into account dependent censoring. Our method extends the factorial analysis of Dobler and Pauly⁷ to the copula-based factorial analysis under dependent censoring. We employ copula models to account for dependent censoring as in Emura and Hsu³¹ who worked on a two-sample test. The proposed method leads to an omnibus test for the global null hypothesis and other tests for composite null hypotheses under factorial designs. The method also allows researchers to conduct sensitivity analyses for treatment effects under a variety of copula models for dependent censoring. We implement the method in an R function surv.factorial(.) in the R package compound.Cox.³²

This article is organized as follows. Section 2 reviews the background for factorial survival analyses. Section 3 proposes a new method for factorial survival analyses under dependent censoring. Section 4 explains the software implementation of the proposed method. Section 5 conducts simulation studies to check the performance of the proposed method. Section 6 provides a real data example. Section 7 concludes the article.

2 Background

We introduce the nonparametric survival analysis method of factorial designs, which was proposed by Dobler and Pauly.⁷ This method defines treatment effects and tests them by using censored survival data collected from one-way, two-way, or other layouts.

2.1 Treatment effects

In clinical trials and animal experiments, the effects of treatments are often assessed by survival outcomes. The classical ANOVA for testing the equality of normal means is inappropriate since survival data are non-normally distributed. Therefore, the nonparametric formulations of treatment effects are needed.^7,33

Consider patients allocated to d different treatments. Let $n_{i}$ be the sample size in the i-th treatment for $i = 1, 2, \dots, d$ . Define the survival function $S_{i} (t) = P (T_{i j} > t)$ for survival time $T_{i j}$ for $j = 1, 2, \dots, n_{i}$ . All survival times are assumed to be independent. Let $N = \sum_{i = 1}^{d} n_{i}$ be the total sample size. The pairwise effect comparing the i-th treatment with the $ℓ$ -th treatment is

\begin{aligned} w_{i ℓ} = P (T_{i 1} > T_{ℓ 1}) + \frac{1}{2} P (T_{i 1} = T_{ℓ 1}) = - \int S_{i}^{\pm} (t) d S_{ℓ} (t) \end{aligned}

(1)

where

S_{i}^{\pm} (t) = {S_{i} (t +) + S_{i} (t -)} / 2

S_{i} (t +) = li m_{s ↓ t} S_{i} (s)

, and

S_{i} (t -) = li m_{s ↑ t} S_{i} (s)

. Note that

w_{i i} = 1 / 2

. This is the probability that a patient in the i-th treatment survives longer than a patient in the

ℓ

-th treatment. The case of

w_{i ℓ} > 1 / 2

(or

w_{i ℓ} < 1 / 2

) gives a beneficial (or harmful) treatment for the i-th treatment relative to the

ℓ

-th treatment. For

d = 2

, this pairwise effects of treatments have been useful in two-sample comparisons.^34–38 However, in factorial analyses involving multiple treatments, the pairwise effects have to be aggregated to characterize the relative treatment effects.

Following the framework of the nonparametric ANOVA,^33,39 Dobler and Pauly⁷ defined the relative treatment effect on the i-th treatment relative to the average treatment effect as

\begin{aligned} p_{i} = \frac{1}{d} \sum_{ℓ = 1}^{d} w_{i ℓ} = - \int S_{i}^{\pm} (t) d \bar{S} (t) \end{aligned}

where

\bar{S} = \sum_{ℓ = 1}^{d} S_{ℓ} / d

. One can write

\begin{aligned} p = (\begin{matrix} p_{1} \\ ⋮ \\ p_{d} \end{matrix}) = A w, w = (w_{11}, \dots, w_{1 d} ⋮ w_{21}, \dots, w_{2 d} ⋮ \dots ⋮ w_{d 1}, \dots, w_{d d})^{'} \end{aligned}

(2)

where

A = I_{d} \otimes ({1^{'}}_{d} / d)

I_{d}

is the

d \times d

identity matrix,

1_{d}^{'} = (1, \dots, 1)

is a

d

-vector of ones, and

\otimes

is the Kronecker product. Hence, the problem of estimating the treatment effects

p

reduces to the problem of estimating the pairwise effects

w

A variety of hypotheses can be formulated using a contrast matrix $C$ such that

\begin{aligned} H_{0} : C p = 0 versus H_{1} : C p \neq 0 \end{aligned}

The simplest example is

C = P_{d} \equiv I_{d} - 1_{d} 1_{d}^{'} / d

, yielding the global hypothesis

\begin{aligned} H_{0} : p_{1} = p_{2} = \dots = p_{d} versus H_{1} : p_{i} \neq p_{ℓ}, \exists (i, ℓ) \end{aligned}

Local hypotheses (e.g.

H_{0} : p_{1} = p_{d}

) can also be formulated by appropriately specifying

C

For a two-way layout, we assume “Factor A” with a levels and “Factor B” with b levels. The treatment effects with interactions are

\begin{aligned} p = (p_{11}, \dots, p_{1 b} ⋮ p_{21}, \dots, p_{2 b} ⋮ \dots ⋮ p_{a 1}, \dots, p_{a b})^{'} \end{aligned}

where

p_{j k}

is the treatment effect for the j-th level in Factor A and the k-th level in Factor B. For instance, the test for main effects in Factor A is

\begin{aligned} H_{0} : p_{1 \cdot} = p_{2 \cdot} = \dots = p_{a \cdot} versus H_{1} : p_{i \cdot} \neq p_{ℓ \cdot} for\; \exists (i, ℓ), i \neq ℓ \end{aligned}

where

p_{i \cdot} = \sum_{ℓ = 1}^{b} p_{i ℓ}

. This hypothesis is specified by

C = P_{a} \otimes 1_{b}^{'} / b

. Other hypotheses can be formulated by choosing appropriate contrast matrices, such as

1_{a}^{'} / a \otimes P_{b}

for main effects in Factor B and

P_{a} \otimes P_{b}

for interactions.

2.2 Factorial analysis under independent censoring

We review the nonparametric test for $H_{0} : C p = 0$ under independent censoring, which was proposed by Dobler and Pauly.⁷ The test rejects $H_{0}$ if $C \hat{p}$ deviates from $0$ , where $\hat{p} = A \hat{w}$ , and $\hat{w}$ is a $d^{2}$ -vector whose element is

\begin{aligned} {\hat{w}}_{i ℓ} = - \int {\hat{S}}_{i}^{\pm} (t) d {\hat{S}}_{ℓ} (t) \end{aligned}

where

{\hat{S}}_{i} (.)

is the KM estimator for

S_{i} (.)

, and

{\hat{S}}_{i}^{\pm} (t) = {{\hat{S}}_{i} (t +) + {\hat{S}}_{i} (t -)} / 2

Let $T = C^{'} (C C^{'})^{+} C$ be a projection matrix, where $({C C}^{'})^{+}$ is the Moore-Penrose inverse of ${C C}^{'}$ that is a unique characterization of the generalized inverse. The F-statistic is defined as

\begin{aligned} F_{N} = \frac{N {\hat{p}}^{'} T \hat{p}}{tr (T \hat{V})} = \frac{N (C \hat{p})^{'} {(C C^{'})}^{+} (C \hat{p})}{tr [{(C C^{'})}^{+} \hat{V} C C^{'}]} \end{aligned}

where

\hat{V}

is the estimator of the asymptotic variance of

\sqrt{N} (\hat{p} - p)

The null distribution of $F_{N}$ was derived by Dobler and Pauly⁷; $F_{N}$ converges to a distribution having the unit mean under $H_{0} : C p = 0$ . Under the alternative hypothesis $H_{1} : C p \neq 0$ , $F_{N}$ goes to infinity in probability. Therefore, a consistent test can be constructed by rejecting $H_{0} : C p = 0$ for $F_{N} > c_{N, α}$ , where $c_{N, α}$ is the upper $α \times 100$ percent point of the null distribution. Dobler and Pauly⁷ suggested calculating the value $c_{N, α}$ via a multiplier bootstrap method. However, the test is valid only under the independent censoring assumption as the KM estimator is inconsistent under dependent censoring.

In this context, our goal is to modify $F_{N}$ to produce a consistent test under dependent censoring.

3. Proposed method under dependent censoring

We first formulate a survival copula model for dependent censoring. Then, we propose a new estimator of treatment effects and a new F-test for treatment effects under dependent censoring.

3.1 Preliminary: Estimating survival under dependent censoring

Recall that survival time $T_{i j}$ may be right-censored by censoring time $U_{i j}$ . The incompletely observable data consist of ${(X_{i j}, δ_{i j}); i = 1, \dots, d, j = 1, \dots, n_{i}}$ , where $X_{i j} = min (T_{i j}, U_{i j})$ and $δ_{i j} = 1 {T_{i j} \leq U_{i j}}$ , where 1{.} is the indicator function.

In order to model the structure of dependent censoring, we postulate the survival copula model:

\begin{aligned} P (T_{i j} > t, U_{i j} > u) = C_{i θ_{ι}} {S_{i} (t), G_{i} (u)}, i = 1, 2, \dots, d \end{aligned}

(3)

where

C_{i θ_{i}} (\cdot, \cdot)

is a copula^40,41 with the parameter

θ_{i}

, and

S_{i} (t) = P (T_{i j} > t)

and

G_{i} (u) = P (U_{i j} > u)

are survival functions whose distributional forms are unspecified. If the Clayton copula⁴² is specified

\begin{aligned} C_{i θ_{i}} (u, v) = (u^{- θ_{i}} + v^{- θ_{i}} - 1)^{- \frac{1}{θ_{i}}}, θ_{i} > 0, i = 1, 2, \dots, d \end{aligned}

where the parameter

θ_{i}

is related to Kendall's tau for

T_{i j}

and

U_{i j}

through

τ_{θ_{i}} = θ_{i} / (θ_{i} + 2)

. If

C_{i θ_{i}} (u, v) = u v

were imposed for

i = 1, 2, \dots, d

, the model would be reduced to the independent censoring model. Other examples of copulas can be seen in Appendix A.1. In many copulas, the limit

θ_{i} \to 0

leads to

C_{i θ_{i}} (u, v) \to u v

. Thus, model (3) provides a quite general dependent censoring regime, including the independent censoring model as its special case of

θ_{i} = 0

. Furthermore, no parametric assumption is required for the distributions of survival time and censoring time.

One can estimate $S_{i} (t)$ using observed data under the assumed copula $C_{i θ_{i}} (\cdot, \cdot)$ in model (3) without any parametric specification for the marginal survival functions. Specifically, we introduce the copula-graphic (CG) estimator,⁴³ in particular, a version of the CG estimator derived by Rivest and Wells³⁰ for the subclass of model (3) given by Archimedean copulas:

\begin{aligned} P (T_{i j} > t, U_{i j} > u) = ϕ_{i θ_{i}}^{- 1} [ϕ_{i θ_{i}} {S_{i} (t)} + ϕ_{i θ_{i}} {G_{i} (u)}], i = 1, \dots, d \end{aligned}

(4)

where

ϕ_{i θ_{i}}

is a generator function that is continuous and strictly decreasing from

ϕ_{i θ_{i}} (0) = \infty

{ϕ_{i}}_{θ_{i}} (1) = 0

.⁴⁰ As for Rivest and Wells,³⁰ we assume that

S_{i} (.)

and

G_{i} (.)

are continuous such that there are no ties in

X_{i j} s

. Then, the CG estimator of Rivest and Wells³⁰ is computed as

\begin{aligned} {\hat{S}}_{i}^{CG} (t) = ϕ_{i θ_{i}}^{- 1} [\sum_{j : X_{i j} \leq t, δ_{i j} = 1} {ϕ_{i θ_{i}} (\frac{{\bar{Y}}_{i} (X_{i j}) - 1}{n_{i}}) - ϕ_{i θ_{i}} (\frac{{\bar{Y}}_{i} (X_{i j})}{n_{i}})}] \end{aligned}

where

{\bar{Y}}_{i} (x) = \sum_{j = 1}^{n_{i}} 1 {X_{i j} \geq x}

is the number at-risk at x. The estimator

{\hat{S}}_{i}^{CG} (.)

is uniformly consistent for

S_{i} (.)

ϕ_{i θ_{i}} (.)

in model (4) is correctly specified: see Theorem 1 of Rivest and Wells.³⁰

Under the Clayton copula, the generator function takes $ϕ_{i θ_{i}} (t) = (t^{- θ_{i}} - 1) / θ_{i}$ for $θ_{i} > 0$ . Then, the CG estimator is

\begin{aligned} {\hat{S}}_{i}^{CG} (t) = {[1 + \sum_{j : X_{i j} \leq t, δ_{i j} = 1} {{(\frac{{\bar{Y}}_{i} (X_{i j}) - 1}{n_{i}})}^{- θ_{i}} - {(\frac{{\bar{Y}}_{i} (X_{i j})}{n_{i}})}^{- θ_{i}}}]}^{- \frac{1}{θ_{i}}}, 0 \leq t \leq max_{j} (X_{i j}) \end{aligned}

Note that the CG estimator reduces to the KM estimator under the independence copula given by

ϕ_{i θ_{i}} (t) = - log (t)

. In this case, one has

\begin{aligned} {\hat{S}}_{i}^{CG} (t) = {\hat{S}}_{i} (t) = \prod_{j : X_{i j} \leq t, δ_{i j} = 1} [1 - \frac{1}{{\bar{Y}}_{i} (X_{i j})}], 0 \leq t \leq max_{j} (X_{i j}) \end{aligned}

The CG estimators under the Gumbel and Frank copulas are given in Appendix A.2.

3.2 Proposed estimator for treatment effects

We shall derive a new estimator for the relative treatment effects $p$ .

We first propose to estimate the pairwise effect $w_{i ℓ}$ in equation (1) by

\begin{aligned} {\hat{w}}_{i ℓ}^{CG} = - \int {\hat{S}}_{i}^{CG \pm} (t) d {\hat{S}}_{ℓ}^{CG} (t) \end{aligned}

Then, by

p = A w

of equation (2), we propose our estimator of

p

\begin{aligned} {\hat{p}}^{CG} = (\begin{matrix} {\hat{p}}_{1}^{CG} \\ ⋮ \\ {\hat{p}}_{d}^{CG} \end{matrix}) = A {\hat{w}}^{CG} = A ({\hat{w}}_{11}^{CG}, \dots, {\hat{w}}_{1 d}^{CG} ⋮ {\hat{w}}_{21}^{CG}, \dots, {\hat{w}}_{2 d}^{CG} ⋮ \dots ⋮ {\hat{w}}_{d 1}^{CG}, \dots, {\hat{w}}_{d d}^{CG})^{'} \end{aligned}

This is a nonparametric estimator in the sense that no parametric assumption is imposed for the marginal distributions.

Under mild conditions and model (4), the estimator ${\hat{p}}^{CG}$ is consistent for $p$ , and permits the asymptotic normal approximation $\sqrt{N} ({\hat{p}}^{CG} - p) \sim N_{d} (0, V),$ $V = A Ω A^{'}$ , where $Ω$ is defined in Appendix A.3. The conditions and proof for the asymptotic normality is given in Appendix A.3.

The variance of ${\hat{p}}^{CG}$ needs to be estimated. However, as the forms of $Ω$ is complex, we suggest applying a jackknife estimator of variance. Following Emura and Hsu,³¹ we estimate the variance of ${\hat{p}}^{CG}$ by

\begin{aligned} V ({\hat{p}}^{CG}) = (\begin{matrix} {\hat{σ}}_{1}^{2} & \dots & {\hat{σ}}_{1 d}^{} \\ ⋮ & ⋱ & ⋮ \\ {\hat{σ}}_{d 1}^{} & \dots & {\hat{σ}}_{d}^{2} \end{matrix}) = \frac{N}{N - 1} \sum_{i = 1}^{d} \sum_{j = 1}^{n_{i}} ({\hat{p}}^{CG (- i j)} - {\hat{p}}^{CG (.)}) ({\hat{p}}^{CG (- i j)} - {\hat{p}}^{CG (.)})^{'} \end{aligned}

where

{\hat{p}}^{CG (- i j)}

is computed without the j-th patient in the i-th treatment, and

{\hat{p}}^{CG (.)} = \sum_{i = 1}^{d} \sum_{j = 1}^{n_{i}} {\hat{p}}^{CG (- i j)} / N

. The standard error (SE) is

SE ({\hat{p}}_{i}^{CG}) = {\hat{σ}}_{i}

, and the

(1 - α) \times 100 %

confidence interval (CI) for

p_{i}

\begin{aligned} [{\hat{p}}_{i}^{CG} - z_{α / 2} SE ({\hat{p}}_{i}^{CG}), {\hat{p}}_{i}^{CG} + z_{α / 2} SE ({\hat{p}}_{i}^{CG})] \end{aligned}

where

z_{α / 2}

is the upper

(α / 2) \times 100 %

point of the standard normal distribution. If necessary, the CI is transformed into [0, 1] by the logit function or others.^7,33,44 Note that the bootstrap resampling method does not work since the CG estimator does not adopt to tied data.

3.3 The proposed test

We shall propose a new test for a linear hypothesis $H_{0} : C p = 0$ versus $H_{1} : C p \neq 0$ under factorial designs, where $C$ is a contrast matrix. The test detects the departure of $C {\hat{p}}^{CG}$ from $0$ under model (4) without any parametric assumption on the marginal survival functions. Let $T = C^{'} (C C^{'})^{+} C$ be a projection matrix. By following the same idea as the F-statistic defined in Section 2.2, we propose a new F-statistic defined as

\begin{aligned} F_{N}^{CG} = \frac{N {\hat{p}}^{CG'} T {\hat{p}}^{CG}}{tr (T {\hat{V}}^{CG})} \end{aligned}

where

{\hat{V}}^{CG} = N \times V ({\hat{p}}^{CG})

is the jackknife estimator of the asymptotic variance of

\sqrt{N} ({\hat{p}}^{CG} - p)

. Then, under

H_{0}

, the distribution of

F_{N}^{CG}

is approximated by a weighted chi-squared distribution

\begin{aligned} F_{N}^{CG} \sim \frac{\sum_{i = 1}^{d} λ_{i} χ_{i}^{2} (1)}{\sum_{i = 1}^{d} λ_{i}} = \frac{\sum_{i = 1}^{d} λ_{i} χ_{i}^{2} (1)}{tr (T V)} \end{aligned}

where

χ_{i}^{2} (f) s

are chi-squared distributed random variables with degrees of freedom f, and

λ_{i} s

are eigenvalues of

T V

. Thus,

F_{N}^{CG}

has the unit mean asymptotically (the proof given in Appendix A.3).

We construct an approximately level $α$ test by rejecting a hypothesis for $F_{N}^{CG} > c_{N, α}$ , where $c_{N, α}$ is the upper $α \times 100$ percentile of the null distribution. Since the critical value $c_{N, α}$ involves unknown quantities, $λ_{i}$ s, we propose two methods to calibrate $c_{N, α}$ :

Simulation method for $c_{N, α}$ : We propose a simulation-based method to calibrate $c_{N, α}$ ; a similar method was suggested by Brunner et al.³³ for data without censoring. This method replaces $λ_{i}$ with ${\hat{λ}}_{i}$ , the eigenvalue of $T {\hat{V}}^{CG}$ . Then, one generates random numbers $χ_{i, (r)}^{2} (1)$ , $i = 1, \dots, d$ , $r = 1, \dots, R$ for the large number R (e.g. $R = 1000$ ). Then, one obtains $c_{N, α}$ as the upper $α \times 100$ percent point of
$\begin{aligned} F_{N, (r)}^{CG} = \frac{\sum_{i = 1}^{d} {\hat{λ}}_{i} χ_{i, (r)}^{2} (1)}{tr (T {\hat{V}}^{CG})}, r = 1, \dots, R \end{aligned}$

Analytical method for $c_{N, α}$ : We propose an analytical method to calibrate $c_{N, α}$ ; a similar method was proposed for data without censoring.^33,39 The idea is to approximate the weighted chi-squared distribution by $g \times χ^{2} (f)$ , where constants g and f are chosen to match the first and second moments.⁴⁵ This leads to
$\begin{aligned} c_{N, α} = \frac{χ^{2} (1 - α, \hat{f})}{\hat{f}}, \hat{f} = \frac{t r^{2} (T {\hat{V}}^{CG})}{tr (T {\hat{V}}^{CG} T {\hat{V}}^{CG})} \end{aligned}$
where $χ^{2} (1 - α, \hat{f})$ is the upper $α$ percentile of the chi-squared distribution with $\hat{f}$ degrees of freedom.

3.4 Restricted follow-up

A technical correction is necessary when survival data are collected within a restricted follow-up time. In this situation, the pairwise effects $w_{i ℓ} = - \int S_{i}^{\pm} (t) d S_{ℓ} (t)$ may not be identifiable from data since the survival function $S_{i} (t)$ is not identifiable for some t and some i. What we can identify from survival data is $S_{i} (t)$ on $t \in [0, τ]$ , where $τ > 0$ is a number satisfying Assumption 2 in Appendix A.3, namely,

\begin{aligned} τ < min_{i} [sup {u : G_{i} (u) S_{i} (u) > 0}] \end{aligned}

As in Dobler and Pauly,⁷ we modify the pairwise effect by

\begin{aligned} w_{i ℓ} & = P [min (T_{i}, τ) > min (T_{ℓ}, τ)] + \frac{1}{2} P [min (T_{i}, τ) = min (T_{ℓ}, τ)] \\ = - \int_{0}^{τ} S_{i}^{\pm} (t) d S_{ℓ} (t) + \frac{S_{i} (τ) S_{ℓ} (τ)}{2} \end{aligned}

The last term on the right-hand side corrects the reduced amount of the pairwise effect due to the truncated integral. Note that

w_{i i} = 1 / 2

. Accordingly, the estimator is also modified as

\begin{aligned} {\hat{w}}_{i ℓ}^{CG} = - \int_{0}^{τ} {\hat{S}}_{i}^{CG \pm} (t) d {\hat{S}}_{ℓ}^{CG} (t) + \frac{{\hat{S}}_{i}^{CG} (τ) {\hat{S}}_{ℓ}^{CG} (τ)}{2} \end{aligned}

3.5 Sensitivity analysis

The proposed method needs to specify a copula (including its parameter $θ_{i}$ ) since it is difficult to be estimated from data.^46–49 Indeed, for the copula to be identifiable from data, one needs to impose proportional hazards models,⁵⁰ parametric marginal survival models,^25,26,29 or other restrictions.⁵¹ If no parametric assumption is made on marginal survival functions, the copula needs to be assumed.^43,52 This is the case for the proposed method.

In this context, we propose a sensitivity analysis under a given copula function with different values of $θ_{i}$ .^{30,31,47–49} If treatment effects are always significant irrespective of the values of $θ_{i}$ , this result confirms the presence of treatment effects under dependent censoring. On the other hand, if the significance depends on the values of $θ_{i}$ , the treatment effects are inconclusive. This is a stringent procedure compared to the traditional one that only examines the treatment effects under the independent censoring model (i.e. $θ_{i} = 0$ ). The usage of this sensitivity analysis is illustrated with a data example in Section 6.

4 Software implementation

We implemented the proposed method in an R function surv.factorial(.) available in the R package compound.Cox.³² The function surv.factorial(.) can compute the estimates of treatment effects ${\hat{p}}^{CG}$ , their SE, and 95% CI (Section 3.2). It also computes the F-statistic and its critical value for testing $H_{0} : C p = 0$ under factorial designs (Section 3.3). The contrast matrix $C$ and the follow-up end $τ$ can be flexibly specified by users. For more details, the readers are referred to the package manual. We have checked the reliability of the R function by extensive simulation studies (Section 5).

For implementation, we assume that the copula is identical across treatments, namely, $C_{1 θ_{1}} = \dots = C_{d θ_{d}}$ and $θ_{1} = \dots = θ_{d}$ . This is because the specifications of d different copulas make the method more difficult to be used and interpreted by users.

5 Simulation

Simulation studies were conducted to investigate the performance of the proposed methods. In particular, we evaluate the estimators ${\hat{p}}^{CG}$ and the test statistic $F_{N}^{CG}$ when the true copula is correctly specified and incorrectly specified (misspecified).

5.1 Simulation designs

The sample size was $n_{i} = 50$ , or 100. We generated $(T_{i j}, U_{i j})$ , $i = 1, 2, \dots, d, j = 1, \dots, n_{i}$ from the Clayton copula model with $θ_{i} = 2$ (Kendall's tau = 0.5) and exponential marginals with the rate parameters $λ_{i}$ for $T_{i j}$ and $μ_{i}$ for $U_{i j}$ . We redefined $U_{i j}$ by $U_{i j} = min (U_{i j}, τ)$ , where $τ = 1$ is the follow-up end. We then obtained the censored data ${(X_{i j}, δ_{i j}); i = 1, \dots, d, j = 1, \dots, n_{i}}$ , where $X_{i j} = min (T_{i j}, U_{i j})$ and $δ_{i j} = 1 {T_{i j} \leq U_{i j}}$ .

Based on the data, we calculated the estimator ${\hat{p}}_{i}^{CG}$ and the 95% CI to see the accuracy of estimating $p_{i}$ . We considered $d = 3$ for the one-way layout (three levels), and $d = 6$ for the two-way layout ( $2 \times 3 = 6$ levels). For the one-way layout, we calculated the F-statistic, $F_{N}^{CG}$ , for testing

\begin{aligned} H_{0} : p_{1} = p_{2} = p_{3} versus H_{1} : p_{1} \neq p_{2}, or p_{1} \neq p_{3}, or p_{2} \neq p_{3} \end{aligned}

by setting the contrast matrix

\begin{aligned} C = I_{3} - \frac{1_{3} 1_{3}^{'}}{3} = \frac{1}{3} (\begin{matrix} 2 & - 1 & - 1 \\ - 1 & 2 & - 1 \\ - 1 & - 1 & 2 \end{matrix}) \end{aligned}

For the two-way layout (

2 \times 3 = 6

levels), the true treatment effects

p

are

\begin{aligned} p = (p_{1}, p_{2}, p_{3}, p_{4}, p_{5}, p_{6})^{'} = {(\underset{First\level\ in\ Factor\A}{\underset{⏟}{p_{11}, p_{12}, p_{13}}}, \underset{Second\level\ in\ Factor\ A}{\underset{⏟}{p_{21}, p_{22}, p_{23}}})}^{'} \end{aligned}

To test the null effect of Factor A, we calculated the F-statistic,

F_{N}^{CG}

, for no main effect in Factor A,

\begin{aligned} H_{0} : p_{11} + p_{12} + p_{13} = p_{21} + p_{22} + p_{23} versus H_{1} : p_{11} + p_{12} + p_{13} \neq p_{21} + p_{22} + p_{23} \end{aligned}

by setting

\begin{aligned} C = P_{2} \otimes \frac{1_{3}^{'}}{3} = \frac{1}{6} (\begin{matrix} \begin{array}{ccc} 1 & 1 & 1 \\ - 1 & - 1 & - 1 \end{array} & \begin{array}{ccc} - 1 & - 1 & - 1 \\ 1 & 1 & 1 \end{array} \end{matrix}) \end{aligned}

The process of estimation and testing did not utilize any knowledge for the exponential marginal distributions (the data-generating process).

The rate parameters $(λ_{i}, μ_{i})$ in the data-generating process were set to yield a realistic amount of treatment effects and censoring percentage $P (U_{i j} < T_{i j}) \times 100$ . We made nine scenarios (Scenarios 1 to 6 for the one-way layout; Scenarios 7 to 9 for the two-way layout) for parameter settings (Table 1). The null hypothesis $H_{0}$ holds for Scenarios 1, 2, 7, and 8 while the alternative hypothesis $H_{1}$ holds for Scenarios 3–6 and 9.

Table 1.
Parameter settings for the data-generating process are defined by nine scenarios. True treatment effects are beneficial (Bold) or harmful (Italic) relative to the overall effect (Normal).

Rate parameters for $T_{i j}$ Rate parameters for $U_{i j}$ True treatment effects

One-way $(λ_{1}, λ_{2}, λ_{3})$ $(μ_{1}, μ_{2}, μ_{3})$ $(p_{1}, p_{2}, p_{3})$

Scenario 1 (1, 1, 1) (1, 1, 1) (0.5, 0.5, 0.5)

Scenario 2 (1, 1, 1) (1, 1.25, 1.5) (0.5, 0.5, 0.5)

Scenario 3 (1, 1.25, 1.5) (1, 1, 1) (0.547, 0.498, 0.455)

Scenario 4 (1, 1.25, 1.5) (1, 1.25, 1.5) (0.547, 0.498, 0.455)

Scenario 5 (1.25, 1, 0.75) (1, 1, 1) (0.447, 0.497, 0.556)

Scenario 6 (1.25, 1, 0.75) (1, 1.25, 1.5) (0.447, 0.497, 0.556)

Two-way $(λ_{11}, λ_{12}, λ_{13}, λ_{21}, λ_{22}, λ_{23})$ $(μ_{11}, μ_{12}, μ_{13}, μ_{21}, μ_{22}, μ_{23})$ $(p_{11}, p_{12}, p_{13}, p_{21}, p_{22}, p_{23})$

Scenario 7 (1, 1, 1, 1, 1, 1) (1, 1, 1, 1, 1, 1) (0.5, 0.5, 0.5, 0.5, 0.5, 0.5)

Scenario 8 (1, 1.25, 1.5, 1, 1.25, 1.5) (1, 1, 1, 1, 1, 1) (0.55,0.50,0.46,0.55,0.50,0.46)

Scenario 9 (1, 1.25, 1.5, 1, 1, 1) (1, 1.25, 1.5, 1, 1, 1) (0.52,0.47,0.43,0.52,0.52,0.52)

	Rate parameters for $T_{i j}$	Rate parameters for $U_{i j}$	True treatment effects
One-way	$(λ_{1}, λ_{2}, λ_{3})$	$(μ_{1}, μ_{2}, μ_{3})$	$(p_{1}, p_{2}, p_{3})$
Scenario 1	(1, 1, 1)	(1, 1, 1)	(0.5, 0.5, 0.5)
Scenario 2	(1, 1, 1)	(1, 1.25, 1.5)	(0.5, 0.5, 0.5)
Scenario 3	(1, 1.25, 1.5)	(1, 1, 1)	(0.547, 0.498, 0.455)
Scenario 4	(1, 1.25, 1.5)	(1, 1.25, 1.5)	(0.547, 0.498, 0.455)
Scenario 5	(1.25, 1, 0.75)	(1, 1, 1)	(0.447, 0.497, 0.556)
Scenario 6	(1.25, 1, 0.75)	(1, 1.25, 1.5)	(0.447, 0.497, 0.556)
Two-way	$(λ_{11}, λ_{12}, λ_{13}, λ_{21}, λ_{22}, λ_{23})$	$(μ_{11}, μ_{12}, μ_{13}, μ_{21}, μ_{22}, μ_{23})$	$(p_{11}, p_{12}, p_{13}, p_{21}, p_{22}, p_{23})$
Scenario 7	(1, 1, 1, 1, 1, 1)	(1, 1, 1, 1, 1, 1)	(0.5, 0.5, 0.5, 0.5, 0.5, 0.5)
Scenario 8	(1, 1.25, 1.5, 1, 1.25, 1.5)	(1, 1, 1, 1, 1, 1)	(0.55,0.50,0.46,0.55,0.50,0.46)
Scenario 9	(1, 1.25, 1.5, 1, 1, 1)	(1, 1.25, 1.5, 1, 1, 1)	(0.52,0.47,0.43,0.52,0.52,0.52)

Based on 1000 replications of simulated data, we assessed the accuracy of the proposed estimators for treatment effects $p$ and proposed tests for $H_{0}$ .

5.2 Simulation results

We first show the results of the estimators and tests when the true value of $θ_{i} = 2$ is known. Next, we turn our attention to the case, where $θ_{i} = 2$ is possibly misspecified as $θ_{i} = 0, 1, 2, 3, and 4$ .

When the true value $θ_{i} = 2$ is known, the means of ${\hat{p}}_{i}^{CG}$ are close to the true value $p_{i}$ , and hence, the estimators are nearly unbiased (Tables A1 and A2 in Supplemental Materials). The variability (standard deviation, SD) of ${\hat{p}}_{i}^{CG}$ vanishes to zero when the sample size gets larger or the censored proportions get lower. This means that the estimators are consistent for the true values. The SDs are very close to the averages of the SEs, which implies the consistency of the jackknife estimates of the asymptotic variance. Consequently, the coverage probabilities of the 95% CI are close to the nominal value of 0.95. Overall, we observe the desirable performance for estimating $p_{i}$ .

Table 2 shows the results for the proposed F-test when the true value of $θ_{i} = 2$ is known. The type I error rates are well-controlled for the nominal level of $α \in {0.10, 0.05, 0.01}$ ; see the rows of “ $H_{0}$ holds”. The proposed test has a reasonable amount of power; see the row of “ $H_{1}$ holds.” The power gets higher when the sample size increases. There is no noticeable difference between the two methods for computing the critical values (simulation method vs. analytical method). Overall, the proposed test has desirable operating characteristics.

Table 2.
Rejection rates for the proposed test based on 1000 simulation runs under level $α$ .

Simulation method^a Analytical method^b

Design Scenario $n_{i}$ $α = 0.10$ $α = 0.05$ $α = 0.01$ $α = 0.10$ $α = 0.05$ $α = 0.01$

$H_{0}$ holds One-way S1 $50$ 0.104 0.056 0.016 0.104 0.056 0.017

$100$ 0.083 0.047 0.008 0.085 0.046 0.008

S2 $50$ 0.102 0.057 0.015 0.102 0.056 0.013

$100$ 0.083 0.040 0.008 0.082 0.040 0.007

Two-way S7 $50$ 0.090 0.053 0.008 0.088 0.055 0.005

$100$ 0.107 0.053 0.016 0.108 0.052 0.017

S8 $50$ 0.092 0.050 0.014 0.091 0.046 0.011

$100$ 0.110 0.051 0.015 0.109 0.052 0.012

$H_{1}$ holds One-way S3 $50$ 0.359 0.246 0.093 0.359 0.250 0.094

$100$ 0.537 0.407 0.192 0.531 0.402 0.179

S4 $50$ 0.329 0.223 0.076 0.334 0.218 0.074

$100$ 0.493 0.360 0.164 0.489 0.358 0.156

S5 $50$ 0.386 0.283 0.131 0.382 0.285 0.129

$100$ 0.652 0.516 0.293 0.646 0.511 0.283

S6 $50$ 0.367 0.252 0.127 0.362 0.253 0.123

$100$ 0.572 0.442 0.228 0.565 0.432 0.231

Two-way S9 $50$ 0.316 0.212 0.077 0.315 0.214 0.072

$100$ 0.493 0.367 0.186 0.492 0.366 0.180

				Simulation method^a	Analytical method^b
$H_{0}$ holds	One-way	S1	$50$	0.104	0.056	0.016	0.104	0.056	0.017
			$100$	0.083	0.047	0.008	0.085	0.046	0.008
		S2	$50$	0.102	0.057	0.015	0.102	0.056	0.013
			$100$	0.083	0.040	0.008	0.082	0.040	0.007
	Two-way	S7	$50$	0.090	0.053	0.008	0.088	0.055	0.005
			$100$	0.107	0.053	0.016	0.108	0.052	0.017
		S8	$50$	0.092	0.050	0.014	0.091	0.046	0.011
			$100$	0.110	0.051	0.015	0.109	0.052	0.012
$H_{1}$ holds	One-way	S3	$50$	0.359	0.246	0.093	0.359	0.250	0.094
			$100$	0.537	0.407	0.192	0.531	0.402	0.179
		S4	$50$	0.329	0.223	0.076	0.334	0.218	0.074
			$100$	0.493	0.360	0.164	0.489	0.358	0.156
		S5	$50$	0.386	0.283	0.131	0.382	0.285	0.129
			$100$	0.652	0.516	0.293	0.646	0.511	0.283
		S6	$50$	0.367	0.252	0.127	0.362	0.253	0.123
			$100$	0.572	0.442	0.228	0.565	0.432	0.231
	Two-way	S9	$50$	0.316	0.212	0.077	0.315	0.214	0.072
			$100$	0.493	0.367	0.186	0.492	0.366	0.180

The critical value $c_{N, α}$ is the upper $α \times 100$ percent point of simulated samples (Section 3.3).

The critical value is $c_{N, α} = χ_{\hat{f}, 1 - α}^{2} / \hat{f}$ with $\hat{f}$ degrees of freedom (Section 3.3).

Figure 1 (the box plots) shows the performance of the estimator ${\hat{p}}_{i}^{CG}$ under the misspecification of $θ_{i}$ with $n_{i} = 100$ . Under Scenario 1, the estimates are nearly unbiased and the coverage probabilities of 95% CI are close to 0.95 across all the misspecified values of $θ_{i}$ . This phenomenon is specific to Scenario 1, where the biases were canceled out by equal censoring percentages across treatments. With unequal censoring percentages as in Scenario 2, the estimates get biased as the specified value of $θ_{i}$ deviates from the true value $θ_{i} = 2$ . Scenario 3 also yields some biases caused by misspecified values of $θ_{i}$ . However, even when $θ_{i}$ is misspecified and the estimator ${\hat{p}}_{i}^{CG}$ is biased, the coverage probabilities remain close to 0.95 in many cases. However, for the misspecified value of $θ_{i} = 0$ , the obvious under-coverage is found.

Figure 1.

Simulation results for ${\hat{p}}_{i}^{CG}, i = 1, 2, 3$ , when $θ_{i} = 2$ (true value) is misspecified as $θ_{i} = 0, 1, 2, 3, 4$ . CP is the coverage probability of the 95% CI. (a) Scenario 1: One-way layout with $p_{1} = p_{2} = p_{3} = 0.500$ . (b) Scenario 2: One-way layout with $p_{1} = p_{2} = p_{3} = 0.500$ . (c) Scenario 3: One-way layout with $p_{1} = 0.547, p_{2} = 0.498, and p_{3} = 0.455$ .

Therefore, if the independent censoring assumption is mistakenly imposed, the estimator ${\hat{p}}_{i}^{CG}$ could be obviously biased and the CI may fail to cover the true value. Similar results are observed under Scenarios 4–6, which are given in Figure B1 of Supplemental Materials. These conclusions for the one-way layout continue to hold for the two-way layout, as seen for Scenarios 7 and 8 (Figures B2 of Supplemental Materials) and Scenario 9 (Figure 2).

Figure 2.

Simulation results for ${\hat{p}}_{i}^{CG}, i = 1, 2, 3$ , when $θ_{i} = 2$ (true value) is misspecified as $θ_{i} = 0, 1, 2, 3, 4$ . A two-way layout is considered under Scenario 9 with $p_{1} = 0.52$ , $p_{2} = 0.47$ , $p_{3} = 0.43$ , and $p_{4} = p_{5} = p_{6} = 0.52$ . CP is the coverage probability of the 95% CI.

6 Data analysis

We analyze breast cancer data available from an R Bioconductor package, curatedBreastData.⁵³ The data contain survival outcomes and other clinical covariates on 2719 breast cancer patients with advanced breast cancer. In the following analysis, we chose a subset of the data having complete information on both disease-free survival (DFS) outcomes and treatment types. There are three treatment types: adjuvant, neoadjuvant, and mixed (both neoadjuvant and adjuvant). Excluding patients with missing DFS or treatment, we shall analyze a subset of 635 patients. We analyze this subset to compare three treatments in the one-way layout ( $n_{1} = 136$ adjuvant, $n_{2} = 107$ mixed, and $n_{3} = 392$ neoadjuvant patients).

Figure 3 displays the estimated survival probabilities for DFS for the three treatments (adjuvant, neoadjuvant, or mixed treatment) by using the CG estimators with the Clayton copula. The adjuvant treatment yields higher survival probabilities than the mixed and neoadjuvant treatments, for all the copula parameters, $θ_{i} = 0$ , $θ_{i} = 2$ , $θ_{i} = 4$ , and $θ_{i} = 8$ . Another obvious pattern is that estimated survival probabilities reduce remarkably when the copula parameter moves from $θ_{i} = 0$ to $θ_{i} = 8$ . This demonstrates the strong impact of dependent censoring on estimated survival.

Figure 3.

Copula-graphic estimates of disease-free survival (DFS) probabilities for the breast cancer data based on three treatments: adjuvant, neoadjuvant, and mixed (both neoadjuvant and adjuvant). The four panels correspond to the fitted Clayton copula with the parameters $θ_{i} = 0$ , $θ_{i} = 2$ , $θ_{i} = 4$ , and $θ_{i} = 8$ .

To check the significance of the difference in treatments, we shall perform a hypothesis test. The log-rank test rejected the global hypothesis $H_{0}^{S} : S_{1} (t) = S_{2} (t) = S_{3} (t) \forall t$ in favor of survival difference (chi-squared statistic = 17.5, critical value = 5.99, P-value < 0.05). We stress that this conclusion was derived under the independent censoring assumption, which corresponds to survival difference at the parameter $θ_{i} = 0$ (see the upper left panel of Figure 3).

In order to see how the conclusion might change under dependent censoring, we applied the proposed F-test. To this end, we perform the test under the parameters $θ_{i} = 0$ , $θ_{i} = 2$ , $θ_{i} = 4$ , and $θ_{i} = 8$ . Table 3 shows that the F-test rejected $H_{0} : p_{1} = p_{2} = p_{3} = 1 / 2$ at a 5% level for all the values of $θ_{i}$ . Hence, even if the independent censoring assumption is violated, the hypothesis $H_{0}^{S} : S_{1} (t) = S_{2} (t) = S_{3} (t) \forall t$ is rejected at a 5% level.

Table 3.

Testing the equality of treatment effects on disease-free survival (DFS) based on the breast cancer data.

		Critical value (5%)	Critical value (5%)
Copula parameter	F-value $F_{N}^{CG}$	$c_{N, α}$ : simulation	$c_{N, α}$ : analytical	P-value
$θ = 0$ (Kendall's tau = 0.00)	133.979	3.073	3.044	<0.0001
$θ = 2$ (Kendall's tau = 0.50)	15.453	3.763	3.515	0.0009
$θ = 4$ (Kendall's tau = 0.67)	22.764	3.602	3.405	0.0001
$θ = 8$ (Kendall's tau = 0.80)	39.807	3.364	3.250	<0.0001

As we detected significant differences in treatment effects, we estimated them by the proposed estimators. Table 4 shows the estimates of $(p_{1}, p_{2}, p_{3})$ using the proposed estimator by fitting the Clayton copula with the parameter $θ_{i}$ . For instance, if $θ_{i} = 8$ is assumed, the adjuvant treatment yields a survival advantage $({\hat{p}}_{1}^{CG} = 0.761 > 0.5)$ over other treatments. The difference is significant since its 95% CI did not cover 0.5 (Table 4). Indeed, for any parameter, the adjuvant treatment yields the best and most significant survival advantage. This strongly confirms the advantage of the adjuvant treatment under dependent censoring. The difference between the neoadjuvant and adjuvant treatments depends on $θ_{i}$ and the comparison is hence inconclusive.

Table 4.

Estimation of treatment effects $p_{i}$ on disease-free survival (DFS) based on the breast cancer data.

Copula parameter	Treatment	Parameter	Estimate	SE	95% CI
$θ = 0$ (Kendall's tau = 0.00)	Adjuvant	$p_{1}$	0.664	0.014	(0.638, 0.691)
	Mixed	$p_{2}$	0.283	0.013	(0.258, 0.307)
	Neoadjuvant	$p_{3}$	0.554	0.015	(0.524, 0.583)
$θ = 2$ (Kendall's tau = 0.50)	Adjuvant	$p_{1}$	0.725	0.034	(0.659, 0.792)
	Mixed	$p_{2}$	0.341	0.032	(0.279, 0.404)
	Neoadjuvant	$p_{3}$	0.436	0.055	(0.329, 0.544)
$θ = 4$ (Kendall's tau = 0.67)	Adjuvant	$p_{1}$	0.754	0.034	(0.686, 0.821)
	Mixed	$p_{2}$	0.364	0.027	(0.310, 0.417)
	Neoadjuvant	$p_{3}$	0.389	0.048	(0.259, 0.482)
$θ = 8$ (Kendall's tau = 0.80)	Adjuvant	$p_{1}$	0.761	0.028	(0.706, 0.815)
	Mixed	$p_{2}$	0.376	0.023	(0.332, 0.421)
	Neoadjuvant	$p_{3}$	0.369	0.035	(0.300, 0.438)

7 Conclusion and discussion

We develop a novel method for factorial survival analysis under a copula-based dependent censoring model. While the majority of the traditional survival analyses were formulated under the independent censoring model, we try to challenge this assumption in a situation of factorial designs. This article is the first attempt to examine dependent censoring in factorial designs. Even though copula-based dependent censoring models have already been considered in a variety of settings and applications,^{12,14,16,17,21,22,27–31} we extend these works toward a new setting: factorial designs.

We derive an F-test that can be applicable to factorial designs with survival outcomes. For the one-way and two-way layouts, the type I error and power properties of the proposed test are well-behaved according to our simulation studies. While we offer two methods to calibrate a critical value of the F-test (the simulation method vs. the analytical method), the practical difference is negligible. The simulation method may have a theoretical advantage due to its more direct approximation to the null distribution, while the analytical method forces a chi-squared distribution to the null distribution. On the other hand, the minor drawback of the simulation method is the need for generating random numbers (i.e. of size $R = 1000$ ) and its relevant random fluctuation of the critical value. However, the required computational time is almost negligible in modern computing environments.

While factorial designs refer to experimentally controlled designs (e.g. randomized clinical trials), the proposed method can also be applicable for observationally collected data with treatment groups as covariates. Indeed, we analyze survival data of 635 breast cancer patients treated by one of three treatments (adjuvant, neoadjuvant, or mixed). The proposed method offers a nonparametric means of estimating treatment effects without assuming any model, such as the proportional hazards model.

Significant survival advantage of the adjuvant treatment is found over the neoadjuvant and mixed treatments. This conclusion is shown to be robust for a variety of dependent censoring scenarios. Nonetheless, the estimated advantage of the adjuvant treatment may be biased since the treatments were not experimentally controlled (i.e. not randomized) and the data were collected observationally. For instance, if the neoadjuvant treatment were administered mainly for patients with poor survival prognosis, their survival would be inherently poor, irrespective of the effect of treatments. Therefore, we cannot exclude the possible bias of the estimated treatment effects due to unbalanced patient characteristics across treatments. It is of interest to develop a method to utilize covariates to adjust for potential biases.

The proposed method relies on the continuity assumption for both survival time and censoring time distributions. Rivest and Wells³⁰ imposed this assumption to derive the CG estimator and its asymptotic properties. For the CG estimator in the proposed method, this assumption is partially relaxed by the restricted follow-up time (Section 3.4). That is, we only need to assume the continuity for $S_{i} (t)$ on $t \in [0, τ]$ . However, the proposed method is not valid for purely discrete survival functions. To see the possible biases due to the violation of the continuity assumption, we conducted additional simulation studies (Supplemental Materials). From the simulation results, we found non-negligible biases in estimation under a discrete survival model. Therefore, it is relevant to extend the proposed method to accommodate discrete survival models.

In our numerical analyses, the Clayton copula was fitted. Nonetheless, other copulas could be tried, such as the Frank copula, especially when negative dependence is suspected for dependent censoring.¹² Our R package (Section 4) allows users to choose the Gumbel and Frank copulas, though we reported our results only for the Clayton copula due to the space limitation. Indeed, it is difficult to justify a suitable copula function by survival data.

The metric of treatment effects is formulated on the basis of pairwise comparison, which is free from model assumptions, such as the proportional hazards model. Recently, pairwise comparison has been adopted for many clinical trials with multiple endpoints⁵⁴ with the aid of Buyse's generalized pairwise comparison.⁵⁵ This approach was extended to survival outcomes with independent censoring by Péron et al.⁵⁶ along with the asymptotic theory.⁵⁷ Extension of the generalized pairwise comparison under dependent censoring is of great interest; see a relevant work for the extension to competing risks outcomes.⁵⁸

Supplemental Material

sj-r-1-smm-10.1177_09622802231215805 - Supplemental material for Factorial survival analysis for treatment effects under dependent censoring

Supplemental material, sj-r-1-smm-10.1177_09622802231215805 for Factorial survival analysis for treatment effects under dependent censoring by Takeshi Emura, Marc Ditzhaus, Dennis Dobler and Kenta Murotani in Statistical Methods in Medical Research

Supplemental Material

sj-r-2-smm-10.1177_09622802231215805 - Supplemental material for Factorial survival analysis for treatment effects under dependent censoring

Supplemental material, sj-r-2-smm-10.1177_09622802231215805 for Factorial survival analysis for treatment effects under dependent censoring by Takeshi Emura, Marc Ditzhaus, Dennis Dobler and Kenta Murotani in Statistical Methods in Medical Research

Supplemental Material

sj-pdf-3-smm-10.1177_09622802231215805 - Supplemental material for Factorial survival analysis for treatment effects under dependent censoring

Supplemental material, sj-pdf-3-smm-10.1177_09622802231215805 for Factorial survival analysis for treatment effects under dependent censoring by Takeshi Emura, Marc Ditzhaus, Dennis Dobler and Kenta Murotani in Statistical Methods in Medical Research

Footnotes

Acknowledgements

The authors thank the Editor and two referees for their valuable suggestions that improved the paper. Emura T was supported financially by JSPS KAKENHI (22K11948; 20H04147). Ditzhaus M was funded by the Deutsche Forschungsgemeinschaft (grant no. DI 2906/1-2). Dobler D would like to thank his new affiliations where a small part of the work has been done: Department of Statistics, TU Dortmund University, and Research Center Trustworthy Data Science and Security, University Alliance Ruhr, Germany. A draft of this manuscript was presented by Emura T in an organized session (EO350: Beyond proportional hazards and standard survival), CM Statistics 2022, London. Comments from the audience helped improve the article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Deutsche Forschungsgemeinschaft, Japan Society for the Promotion of Science (grant number DI 2906/1-2, 20H04147, 22K11948).

ORCID iDs

Takeshi Emura

Dennis Dobler

Supplementary data

Online supplemental materials related to this article are additional simulation results, the R code for the simulation studies, and the R code for the data analysis.

Appendices

References

Montgomery

. Design and analysis of experiments. 10th ed. Hoboken, NJ: John Wiley & Sons, 2019.

Mishra

Singh

Pandey

, et al. Application of student's t-test, analysis of variance, and covariance. Ann Card Anaesth 2019; 22: 407.

Futschik

Taus

Zehetmayer

. An omnibus test for the global null hypothesis. Stat Methods Med Res 2019; 28: 2292–2304.

Fizazi

Foulon

Carles

, et al. Abiraterone plus prednisone added to androgen deprivation therapy and docetaxel in de novo metastatic castration-sensitive prostate cancer (PEACE-1): a multicentre, open-label, randomised, phase 3 study with a 2 × 2 factorial design. Lancet 2022; 399: 1695–1707.

Freidlin

Korn

. Two-by-two factorial cancer treatment trials: is sufficient attention being paid to possible interactions? J Natl Cancer Inst 2017; 109: djx146.

Tewari

Sill

, et al. Improved survival with bevacizumab in advanced cervical cancer. N Engl J Med 2014; 370: 734–743.

Dobler

Pauly

. Factorial analyses of treatment effects under independent right-censoring. Stat Methods Med Res 2020; 29: 325–343.

Gorfine

Schlesinger

Hsu

. K-sample omnibus non-proportional hazards tests based on right-censored data. Stat Methods Med Res 2020; 29: 2830–2850.

Ditzhaus

Dobler

Pauly

. Inferring median survival differences in general factorial designs via permutation tests. Stat Methods Med Res 2021; 30: 875–891.

10.

Ditzhaus

Genuneit

Janssen

, et al. CASANOVA: permutation inference in factorial survival designs. Biometrics 2023; 79: 203–215.

11.

Klein

Moeschberger

. Survival analysis: techniques for censored and truncated data. New York: Springer, 2003.

12.

Braekers

Veraverbeke

. A copula-graphic estimator for the conditional survival function under dependent censoring. Can J Stat 2005; 33: 429–447.

13.

Zhang

. Semiparametric likelihood estimation in survival models with informative censoring. J Multivar Anal 2012; 106: 187–211.

14.

Emura

Michimae

. A copula-based inference to piecewise exponential models under dependent censoring, with application to time to metamorphosis of salamander larvae. Environ Ecol Stat 2017; 24: 151–173.

15.

Staplin

Kimber

Collett

, et al. Dependent censoring in piecewise exponential survival models. Stat Methods Med Res 2015; 24: 325–341.

16.

Emura

Chen

. Gene selection for survival data under dependent censoring, a copula-based approach. Stat Methods Med Res 2016; 25: 2840–2857.

17.

Emura

Chen

. Analysis of survival data with dependent censoring: copula-based approaches. Singapore: Springer, 2018.

18.

Dettoni

Marra

Radice

. Generalized link-based additive survival models with informative censoring. J Comput Graph Stat 2020; 29: 503–512.

19.

Andersen

Perme

. Pseudo-observations in survival analysis. Stat Methods Med Res 2010; 19: 71–99.

20.

Collett

. Modelling survival data in medical research fourth edition. Oxfordshire: Chapman & Hall, 2023.

21.

Connors

, et al. Proportional hazard model estimation under dependent censoring using copulas and penalized likelihood. Stat Med 2018; 37: 2238–2251.

22.

Moradian

Larocque

Bellavance

. Survival forests for data with dependent censoring. Stat Methods Med Res 2019; 28: 445–461.

23.

Schneider

Demarqui

Colosimo

, et al. An approach to model clustered survival data with dependent censoring. Biom J 2020; 62: 157–174.

24.

Schneider

Demarqui

de Freitas Costa

. Free-ranging dogs’ lifetime estimated by an approach for long-term survival data with dependent censoring. Environ Ecol Stat 2022; 29: 869–911.

25.

Deresa

Van Keilegom

. Flexible parametric model for survival data subject to dependent censoring. Biom J 2020; 62: 136–156.

26.

Deresa

Van Keilegom

. A multivariate normal regression model for survival data subject to different types of dependent censoring. Comput Stat Data Anal 2020; 144: 106879.

27.

Wang

. Evaluating association between two event times with observations subject to informative censoring. J Am Stat Assoc 2021; 118: 542.

28.

Deresa

Van Keilegom

Antonio

. Copula-based inference for bivariate survival data with left truncation and dependent censoring. Insur Math Econ 2022; 107: 1–21.

29.

Czado

Van Keilegom

. Dependent censoring based on parametric copulas. Biometrika 2023; 110: 721–738.

30.

Rivest

Wells

. A martingale approach to the copula-graphic estimator for the survival function under dependent censoring. J Multivar Anal 2001; 79: 138–155.

31.

Emura

Hsu

. Estimation of the Mann–Whitney effect in the two-sample problem under dependent censoring. Comput Stat Data Anal 2020; 150: 106990.

32.

Emura

Matsui

Chen

. Compound.cox: univariate feature selection and compound covariate for predicting survival. Comput Methods Programs Biomed 2019; 168: 21–37.

33.

Brunner

Konietschke

Pauly

, et al. Rank-based procedures in factorial designs: hypotheses about non-parametric treatment effects. J R Stat Soc: Ser B 2017; 79: 1463–1485.

34.

Efron

. The two sample problem with censored data. In: Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, Vol. 4, 1967, pp.831–853. Berkeley, CA: University of California Press.

35.

Koziol

Jia

. The concordance index C and the Mann–Whitney parameter Pr (X>Y) with randomly censored data. Biom J 2009; 51: 467–474.

36.

Dobler

Pauly

. Bootstrap- and permutation-based inference for the Mann–Whitney effect for right-censored and tied data. Test 2018; 27: 639–658.

37.

Dobler

Friedrich

Pauly

. Nonparametric MANOVA in meaningful effects. Ann Inst Stat Math 2020; 72: 997–1022.

38.

Nowak

Mütze

Konietschke

. Group sequential methods for the Mann-Whitney parameter. Stat Methods Med Res 2022; 31: 2004–2020.

39.

Brunner

Puri

. Nonparametric methods in factorial designs. Stat Pap 2001; 42: 1–52.

40.

Nelsen

. An introduction to copulas. New York, NY: Springer Science & Business Media, 2006.

41.

Durante

Sempi

. Principles of copula theory. Boca Raton, FL: CRC Press, 2016.

42.

Clayton

. A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence. Biometrika 1978; 65: 141–151.

43.

Zheng

Klein

. Estimates of marginal survival for dependent competing risks based on an assumed copula. Biometrika 1995; 82: 127–138.

44.

Perme

Manevski

. Confidence intervals for the Mann–Whitney test. Stat Methods Med Res 2019; 28: 3755–3768.

45.

Box

. Some theorems on quadratic forms applied in the study of analysis of variance problems, I. Effect of inequality of variance in the one-way classification. Ann Math Stat 1954; 25: 290–302.

46.

Tsiatis

. A nonidentifiability aspect of the problem of competing risks. Proc Natl Acad Sci USA 1975; 72: 20–22.

47.

De Uña-Álvarez

Veraverbeke

. Generalized copula-graphic estimator. Test 2013; 22: 343–360.

48.

De Uña-Álvarez

Veraverbeke

. Copula-graphic estimation with left-truncated and right-censored data. Statistics (Ber) 2017; 51: 387–403.

49.

Yeh

Liao

Emura

. Sensitivity analysis for survival prognostic prediction with gene selection: a copula method for dependent censoring. Biomedicines 2023; 11: 797.

50.

Heckman

Honoré

. The identifiability of the competing risks model. Biometrika 1989; 76: 325–330.

51.

Schwarz

Jongbloed

Van Keilegom

. On the identifiability of copulas in bivariate competing risks models. Can J Stat 2013; 41: 291–303.

52.

Carrière

. Removing cancer when it is correlated with other causes of death. Biom J 1995; 37: 339–350.

53.

Planey

. curatedBreastData: curated breast cancer gene expression data with survival and treatment information. CRAN , R package version 2.26.0, 2022.

54.

Deltuvaite-Thomas

Verbeeck

Burzykowski

, et al. Generalized pairwise comparisons for censored data: an overview. Biom J 2023; 65: 2100354.

55.

Buyse

. Generalized pairwise comparisons of prioritized outcomes in the two-sample problem. Stat Med 2010; 29: 3245–3257.

56.

Péron

Buyse

Ozenne

, et al. An extension of generalized pairwise comparisons for prioritized outcomes in the presence of censoring. Stat Methods Med Res 2018; 27: 1230–1239.

57.

Ozenne

Budtz-Jørgensen

Péron

. The asymptotic distribution of the net benefit estimator in presence of right-censoring. Stat Methods Med Res 2021; 30: 2399–2412.

58.

Cantagallo

De Backer

Kicinski

, et al. A new measure of treatment effect in clinical trials involving competing risks based on generalized pairwise comparisons. Biom J 2021; 63: 272–288.

59.

Gumbel

. Distributions des valeurs extremes en plusiers dimensions. Publ Inst Stat Univ Paris 1960; 9: 171–173.

60.

Peng

Xiang

Wang

. Semiparametric regression analysis of clustered survival data with semi-competing risks. Comput Stat Data Anal 2018; 124: 53–70.

61.

Frank

. On the simultaneous associativity of F(x,y) and x + y – F(x,y). Aequationes Math 1979; 19: 194–226.

62.

de Oliveira Peres

Achcar

Martinez

. Bivariate lifetime models in presence of cure fraction: a comparative study with many different copula functions. Heliyon 2020; 6: e03961.

63.

Huang

Wang

Emura

. A copula-based Markov chain model for serially dependent event times with a dependent terminal event. Jpn J Stat Data Sci 2021; 4: 917–951.

64.

Morgenstern

. Einfache beispiele zweidimensionaler verteilungen. Mitteilingsblatt Math Stat 1956; 8: 234–235.

65.

Ota

Kimura

. Effective estimation algorithm for parameters of multivariate Farlie–Gumbel–Morgenstern copula. Jpn J Stat Data Sci 2021; 4: 1049–1078.

66.

Susam

. A multi-parameter generalized Farlie-Gumbel-Morgenstern bivariate copula family via Bernstein polynomial. Hacet J Math Stat 2022; 51: 618–631.

67.

de Oliveira

de Oliveira Peres

dos Santos

, et al. A Bayesian inference approach for bivariate Weibull distributions derived from Roy and Morgenstern methods. Stat Optim Inf Comput 2021; 9: 529–554.

68.

Wilke

. A copula model for dependent competing risks. J R Stat Soc: Ser C (Appl Stat) 2010; 59: 359–376.

69.

Ditzhaus

Smaga

. Permutation test for the multivariate coefficient of variation in factorial designs. J Multivar Anal 2022; 187: 104848.

70.

Billingsley

. Convergence of probability measures. New York, NY: John Wiley & Sons, 1999.

71.

Van der Vaart

. Asymptotic statistics. Cambridge: Cambridge University Press, 2000.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.01 MB

1.33 MB

				Simulation method^a			Analytical method^b
	Design	Scenario	$n_{i}$	$α = 0.10$	$α = 0.05$	$α = 0.01$	$α = 0.10$	$α = 0.05$	$α = 0.01$
$H_{0}$ holds	One-way	S1	$50$	0.104	0.056	0.016	0.104	0.056	0.017
			$100$	0.083	0.047	0.008	0.085	0.046	0.008
		S2	$50$	0.102	0.057	0.015	0.102	0.056	0.013
			$100$	0.083	0.040	0.008	0.082	0.040	0.007
	Two-way	S7	$50$	0.090	0.053	0.008	0.088	0.055	0.005
			$100$	0.107	0.053	0.016	0.108	0.052	0.017
		S8	$50$	0.092	0.050	0.014	0.091	0.046	0.011
			$100$	0.110	0.051	0.015	0.109	0.052	0.012
$H_{1}$ holds	One-way	S3	$50$	0.359	0.246	0.093	0.359	0.250	0.094
			$100$	0.537	0.407	0.192	0.531	0.402	0.179
		S4	$50$	0.329	0.223	0.076	0.334	0.218	0.074
			$100$	0.493	0.360	0.164	0.489	0.358	0.156
		S5	$50$	0.386	0.283	0.131	0.382	0.285	0.129
			$100$	0.652	0.516	0.293	0.646	0.511	0.283
		S6	$50$	0.367	0.252	0.127	0.362	0.253	0.123
			$100$	0.572	0.442	0.228	0.565	0.432	0.231
	Two-way	S9	$50$	0.316	0.212	0.077	0.315	0.214	0.072
			$100$	0.493	0.367	0.186	0.492	0.366	0.180