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.

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_{ij}>t)$ for survival time $T_{ij}$ 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 $\ell$ -th treatment is

$$\begin{array}{r}{w}_{i\ell }=P(T_{i1}>T_{\ell 1})+\frac{1}{2}P(T_{i1}=T_{\ell 1})=-\int S_i^{\pm }(t)d{S}_{\ell }(t)\end{array}$$

(1)

where $S_i^{\pm }(t)=\{S_i(t+)+S_i(t-)\}/2$ , $S_i(t+)=\text{li}{\text{m}}_{s\downarrow t}{S}_i(s)$ , and $S_i(t-)=\text{li}{\text{m}}_{s\uparrow t}{S}_i(s)$ . Note that $w_{ii}=1/2$ . This is the probability that a patient in the i-th treatment survives longer than a patient in the $\ell$ -th treatment. The case of $w_{i\ell }>1/2$ (or $w_{i\ell }<1/2$ ) gives a beneficial (or harmful) treatment for the i-th treatment relative to the $\ell$ -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{array}{r}{p}_i=\frac{1}{d}\sum _{\ell =1}^d{w}_{i\ell }=-\int S_i^{\pm }(t)d\overline{S}(t)\end{array}$$

where $\overline{S}={\sum }_{\ell =1}^d{S}_{\ell }/d$ . One can write

$$\begin{array}{r}\mathit{p}=\left(\begin{array}{c}{p}_1\\ ⋮\\ p_d\end{array}\right)=\mathit{A}\mathit{w},\mathit{w}=(w_{11},\dots ,w_{1d}⋮w_{21},\dots ,w_{2d}⋮\dots ⋮w_{d1},\dots ,w_{dd}{)}^{\mathrm{\prime }}\end{array}$$

(2)

where $\mathit{A}={\mathit{I}}_d\otimes ({{\mathbf{\text{1}}}^{\mathrm{\prime }}}_d/d)$ , ${\mathit{I}}_d$ is the $d\times d$ identity matrix, ${\mathbf{\text{1}}}_d^{\mathrm{\prime }}=(1,\dots ,1)$ is a $d$ -vector of ones, and $\otimes$ is the Kronecker product. Hence, the problem of estimating the treatment effects $\mathit{p}$ reduces to the problem of estimating the pairwise effects  $\mathit{w}$ .

A variety of hypotheses can be formulated using a contrast matrix $\mathit{C}$ such that

$$\begin{array}{r}{H}_0:\mathit{C}\mathit{p}=\mathbf{\text{0}}\text{versus}{H}_1:\mathit{C}\mathit{p}\ne \mathbf{\text{0}}\end{array}$$

The simplest example is $\mathit{C}={\mathit{P}}_d\equiv {\mathit{I}}_d-{\mathbf{\text{1}}}_d{\mathbf{\text{1}}}_d^{\mathrm{\prime }}/\text{d}$ , yielding the global hypothesis

$$\begin{array}{r}{H}_0:p_1=p_2=\cdots =p_d\text{versus}{H}_1:p_i\ne p_{\ell },\mathrm{\exists }(i,\ell )\end{array}$$

Local hypotheses (e.g. $H_0:p_1=p_d$ ) can also be formulated by appropriately specifying $\mathit{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{array}{r}\mathit{p}=(p_{11},\dots ,p_{1b}⋮p_{21},\dots ,p_{2b}⋮\dots ⋮p_{a1},\dots ,p_{ab}{)}^{\mathrm{\prime }}\end{array}$$

where $p_{jk}$ 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{array}{r}{H}_0:p_{1\cdot }=p_{2\cdot }=\cdots =p_{a\cdot }\text{versus}{H}_1:p_{i\cdot }\ne p_{\ell \cdot }\text{for; }\mathrm{\exists }(i,\ell ),i\ne \ell \end{array}$$

where $p_{i\cdot }={\sum }_{\ell =1}^b{p}_{i\ell }$ . This hypothesis is specified by $\mathit{C}={\mathit{P}}_a\otimes {\mathbf{\text{1}}}_b^{\mathrm{\prime }}/b$ . Other hypotheses can be formulated by choosing appropriate contrast matrices, such as ${\mathbf{\text{1}}}_a^{\mathrm{\prime }}/a\otimes {\mathit{P}}_b$ for main effects in Factor B and ${\mathit{P}}_a\otimes {\mathit{P}}_b$ for interactions.

2.2 Factorial analysis under independent censoring

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

$$\begin{array}{r}{\hat{w}}_{i\ell }=-\int {\hat{S}}_i^{\pm }(t)d{\hat{S}}_{\ell }(t)\end{array}$$

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 $\mathit{T}={\mathit{C}}^{\mathbf{\prime }}(\mathit{C}{\mathit{C}}^{\mathbf{\prime }}{)}^{+}\mathit{C}$ be a projection matrix, where $({\mathit{C}\mathit{C}}^{\mathrm{\prime }}{)}^{+}$ is the Moore-Penrose inverse of ${\mathit{C}\mathit{C}}^{\mathrm{\prime }}$ that is a unique characterization of the generalized inverse. The F-statistic is defined as

$$\begin{array}{r}{F}_N=\frac{N{\hat{\mathit{p}}}^{\mathrm{\prime }}\mathit{T}\hat{\mathit{p}}}{\text{tr}(\mathit{T}\hat{\mathit{V}})}=\frac{N(\mathit{C}\hat{\mathit{p}}){}^{\mathbf{\prime }}{(\mathit{C}{\mathit{C}}^{\mathbf{\prime }})}^{+}(\mathit{C}\hat{\mathit{p}})}{\text{tr}[{(\mathit{C}{\mathit{C}}^{\mathbf{\prime }})}^{+}\hat{\mathit{V}}\mathit{C}{\mathit{C}}^{\mathbf{\prime }}]}\end{array}$$

where $\hat{\mathit{V}}$ is the estimator of the asymptotic variance of $\sqrt{N}(\hat{\mathit{p}}-\mathit{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:\mathit{C}\mathit{p}=\mathbf{\text{0}}$ . Under the alternative hypothesis $H_1:\mathit{C}\mathit{p}\ne \mathbf{\text{0}}$ , $F_N$ goes to infinity in probability. Therefore, a consistent test can be constructed by rejecting $H_0:\mathit{C}\mathit{p}=\mathbf{\text{0}}$ for $F_N>c_{N,\alpha }$ , where $c_{N,\alpha }$ is the upper $\alpha \times 100$ percent point of the null distribution. Dobler and Pauly ⁷ suggested calculating the value $c_{N,\alpha }$ 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_{ij}$ may be right-censored by censoring time $U_{ij}$ . The incompletely observable data consist of $\{(X_{ij},{\delta }_{ij});i=1,\dots ,d,j=1,\dots ,n_i\}$ , where $X_{ij}=\text{min}(T_{ij},U_{ij})$ and ${\delta }_{ij}=1\{T_{ij}\le U_{ij}\}$ , where 1{.} is the indicator function.

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

$$\begin{array}{r}P(T_{ij}>t,U_{ij}>u)=C_{i{\theta }_{\iota }}\{S_i(t),G_i(u)\},i=1,2,\dots ,d\end{array}$$

(3)

where $C_{i{\theta }_i}(\cdot ,\cdot )$ is a copula^40,41 with the parameter ${\theta }_i$ , and $S_i(t)=P(T_{ij}>t)$ and $G_i(u)=P(U_{ij}>u)$ are survival functions whose distributional forms are unspecified. If the Clayton copula ⁴² is specified

$$\begin{array}{r}{C}_{i{\theta }_i}(u,v)=(u^{-{\theta }_i}+v^{-{\theta }_i}-1{)}^{-\frac{1}{{\theta }_i}},{\theta }_i>0,i=1,2,\dots ,d\end{array}$$

where the parameter ${\theta }_i$ is related to Kendall's tau for $T_{ij}$ and $U_{ij}$ through ${\tau }_{{\theta }_i}={\theta }_i/({\theta }_i+2)$ . If $C_{i{\theta }_i}(u,v)=uv$ 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 ${\theta }_i\to 0$ leads to $C_{i{\theta }_i}(u,v)\to uv$ . Thus, model (3) provides a quite general dependent censoring regime, including the independent censoring model as its special case of ${\theta }_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{\theta }_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{array}{r}P(T_{ij}>t,U_{ij}>u)={\varphi }_{i{\theta }_i}^{-1}[{\varphi }_{i{\theta }_i}\{S_i(t)\}+{\varphi }_{i{\theta }_i}\{G_i(u)\}],i=1,\dots ,d\end{array}$$

(4)

where ${\varphi }_{i{\theta }_i}$ is a generator function that is continuous and strictly decreasing from ${\varphi }_{i{\theta }_i}(0)=\mathrm{\infty }$ to ${{\varphi }_i}_{{\theta }_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_{ij}\text{s}$ . Then, the CG estimator of Rivest and Wells ³⁰ is computed as

$$\begin{array}{r}{\hat{S}}_i^{\text{CG}}(t)={\varphi }_{i{\theta }_i}^{-1}\left[\sum _{j:X_{ij}\le t,{\delta }_{ij}=1}\left\{{\varphi }_{i{\theta }_i}\left(\frac{{\overline{Y}}_i(X_{ij})-1}{n_i}\right)-{\varphi }_{i{\theta }_i}\left(\frac{{\overline{Y}}_i(X_{ij})}{n_i}\right)\right\}\right]\end{array}$$

where ${\overline{Y}}_i(x)={\sum }_{j=1}^{n_i}1\{X_{ij}\ge x\}$ is the number at-risk at x. The estimator ${\hat{S}}_i^{\text{CG}}(.)$ is uniformly consistent for $S_i(.)$ if ${\varphi }_{i{\theta }_i}(.)$ in model (4) is correctly specified: see Theorem 1 of Rivest and Wells. ³⁰

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

$$\begin{array}{r}{\hat{S}}_i^{\text{CG}}(t)={\left[1+\sum _{j:X_{ij}\le t,{\delta }_{ij}=1}\left\{{\left(\frac{{\overline{Y}}_i(X_{ij})-1}{n_i}\right)}^{-{\theta }_i}-{\left(\frac{{\overline{Y}}_i(X_{ij})}{n_i}\right)}^{-{\theta }_i}\right\}\right]}^{-\frac{1}{{\theta }_i}},0\le t\le \underset{j}{\text{max}}(X_{ij})\end{array}$$

Note that the CG estimator reduces to the KM estimator under the independence copula given by ${\varphi }_{i{\theta }_i}(t)=-\text{log}(t)$ . In this case, one has

$$\begin{array}{r}{\hat{S}}_i^{\text{CG}}(t)={\hat{S}}_i(t)=\prod _{j:X_{ij}\le t,{\delta }_{ij}=1}\left[1-\frac{1}{{\overline{Y}}_i(X_{ij})}\right],0\le t\le \underset{j}{\text{max}}(X_{ij})\end{array}$$

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 $\mathit{p}$ .

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

$$\begin{array}{r}{\hat{w}}_{i\ell }^{\text{CG}}=-\int {\hat{S}}_i^{\text{CG}\pm }(t)d{\hat{S}}_{\ell }^{\text{CG}}(t)\end{array}$$

Then, by $\mathit{p}=\mathit{A}\mathit{w}$ of equation (2), we propose our estimator of $\mathit{p}$ as

$$\begin{array}{r}{\hat{\mathit{p}}}^{\text{CG}}=\left(\begin{array}{c}{\hat{p}}_1^{\text{CG}}\\ ⋮\\ {\hat{p}}_d^{\text{CG}}\end{array}\right)=\mathit{A}{\hat{\mathit{w}}}^{\text{CG}}=\mathit{A}({\hat{w}}_{11}^{\text{CG}},\dots ,{\hat{w}}_{1d}^{\text{CG}}⋮{\hat{w}}_{21}^{\text{CG}},\dots ,{\hat{w}}_{2d}^{\text{CG}}⋮\dots ⋮{\hat{w}}_{d1}^{\text{CG}},\dots ,{\hat{w}}_{dd}^{\text{CG}}{)}^{\mathrm{\prime }}\end{array}$$

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{\mathit{p}}}^{\text{CG}}$ is consistent for $\mathit{p}$ , and permits the asymptotic normal approximation $\sqrt{N}({\hat{\mathit{p}}}^{\text{CG}}-\mathit{p})\sim N_d(0,\mathit{V}),$ $\mathit{V}=\mathit{A}\mathbf{\Omega }{\mathit{A}}^{\mathbf{\prime }}$ , where $\mathbf{\Omega }$ is defined in Appendix A.3. The conditions and proof for the asymptotic normality is given in Appendix A.3.

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

$$\begin{array}{r}\mathit{V}({\hat{\mathit{p}}}^{\text{CG}})=\left(\begin{array}{ccc}{\hat{\sigma }}_1^2& \dots & {\hat{\sigma }}_{1d}^{}\\ ⋮& \ddots & ⋮\\ {\hat{\sigma }}_{d1}^{}& \dots & {\hat{\sigma }}_d^2\end{array}\right)=\frac{N}{N-1}\sum _{i=1}^d\sum _{j=1}^{n_i}({\hat{\mathit{p}}}^{\text{CG}(-ij)}-{\hat{\mathit{p}}}^{\text{CG}(.)})({\hat{\mathit{p}}}^{\text{CG}(-ij)}-{\hat{\mathit{p}}}^{\text{CG}(.)}){}^{\mathbf{\prime }}\end{array}$$

where ${\hat{\mathit{p}}}^{\text{CG}(-ij)}$ is computed without the j-th patient in the i-th treatment, and ${\hat{\mathit{p}}}^{\text{CG}(.)}={\sum }_{i=1}^d{\sum }_{j=1}^{n_i}{\hat{\mathit{p}}}^{\text{CG}(-ij)}/N$ . The standard error (SE) is $\text{SE}({\hat{p}}_i^{\text{CG}})={\hat{\sigma }}_i$ , and the $(1-\alpha )\times 100\mathrm{\%}$ confidence interval (CI) for $p_i$ is

$$\begin{array}{r}[{\hat{p}}_i^{\text{CG}}-z_{\alpha /2}\text{SE}({\hat{p}}_i^{\text{CG}}),{\hat{p}}_i^{\text{CG}}+z_{\alpha /2}\text{SE}({\hat{p}}_i^{\text{CG}})]\end{array}$$

where $z_{\alpha /2}$ is the upper $(\alpha /2)\times 100\mathrm{\%}$ 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:\mathit{C}\mathit{p}=\mathbf{\text{0}}$ versus $H_1:\mathit{C}\mathit{p}\ne \mathbf{\text{0}}$ under factorial designs, where $\mathit{C}$ is a contrast matrix. The test detects the departure of $\mathit{C}{\hat{\mathit{p}}}^{\text{CG}}$ from $\mathbf{\text{0}}$ under model (4) without any parametric assumption on the marginal survival functions. Let $\mathit{T}={\mathit{C}}^{\mathbf{\prime }}(\mathit{C}{\mathit{C}}^{\mathbf{\prime }}{)}^{+}\mathit{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{array}{r}{F}_N^{\text{CG}}=\frac{N{\hat{\mathit{p}}}^{\text{CG}\mathrm{\prime }}\mathit{T}{\hat{\mathit{p}}}^{\text{CG}}}{\text{tr}(\mathit{T}{\hat{\mathit{V}}}^{\text{CG}})}\end{array}$$

where ${\hat{\mathit{V}}}^{\text{CG}}=N\times \mathit{V}({\hat{\mathit{p}}}^{\text{CG}})$ is the jackknife estimator of the asymptotic variance of $\sqrt{N}({\hat{\mathit{p}}}^{\text{CG}}-\mathit{p})$ . Then, under $H_0$ , the distribution of $F_N^{\text{CG}}$ is approximated by a weighted chi-squared distribution

$$\begin{array}{r}{F}_N^{\text{CG}}\sim \frac{{\sum }_{i=1}^d{\lambda }_i{\chi }_i^2(1)}{{\sum }_{i=1}^d{\lambda }_i}=\frac{{\sum }_{i=1}^d{\lambda }_i{\chi }_i^2(1)}{\text{tr}(\mathit{T}\mathit{V})}\end{array}$$

where ${\chi }_i^2(f)\text{s}$ are chi-squared distributed random variables with degrees of freedom f, and ${\lambda }_i\text{s}$ are eigenvalues of $\mathit{T}\mathit{V}$ . Thus, $F_{\text{N}}^{\text{CG}}$ has the unit mean asymptotically (the proof given in Appendix A.3).

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

Simulation method for $c_{N,\alpha }$ : We propose a simulation-based method to calibrate $c_{N,\alpha }$ ; a similar method was suggested by Brunner et al. ³³ for data without censoring. This method replaces ${\lambda }_i$ with ${\hat{\lambda }}_i$ , the eigenvalue of $\mathit{T}{\hat{\mathit{V}}}^{\text{CG}}$ . Then, one generates random numbers ${\chi }_{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,\alpha }$ as the upper $\alpha \times 100$ percent point of

$$\begin{array}{r}{F}_{N,(r)}^{\text{CG}}=\frac{{\sum }_{i=1}^d{\hat{\lambda }}_i{\chi }_{i,(r)}^2(1)}{\text{tr}(\mathit{T}{\hat{\mathit{V}}}^{\text{CG}})},r=1,\dots ,R\end{array}$$

Analytical method for $c_{N,\alpha }$ : We propose an analytical method to calibrate $c_{N,\alpha }$ ; a similar method was proposed for data without censoring.^33,39 The idea is to approximate the weighted chi-squared distribution by $g\times {\chi }^2(f)$ , where constants g and f are chosen to match the first and second moments. ⁴⁵ This leads to

$$\begin{array}{r}{c}_{N,\alpha }=\frac{{\chi }^2(1-\alpha ,\hat{f})}{\hat{f}},\hat{f}=\frac{\text{t}{\text{r}}^2(\mathit{T}{\hat{\mathit{V}}}^{\text{CG}})}{\text{tr}(\mathit{T}{\hat{\mathit{V}}}^{\text{CG}}\mathit{T}{\hat{\mathit{V}}}^{\text{CG}})}\end{array}$$

where ${\chi }^2(1-\alpha ,\hat{f})$ is the upper $\alpha$ 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\ell }=-\int S_i^{\pm }(t)d{S}_{\ell }(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,\tau ]$ , where $\tau >0$ is a number satisfying Assumption 2 in Appendix A.3, namely,

$$\begin{array}{r}\tau <\underset{i}{min}[sup\{u:G_i(u)S_i(u)>0\}]\end{array}$$

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

$$\begin{array}{rl}{w}_{i\ell }& =P[\text{min}(T_i,\tau )>\text{min}(T_{\ell },\tau )]+\frac{1}{2}P[\text{min}(T_i,\tau )=\text{min}(T_{\ell },\tau )]\\ & =-\underset{0}{\overset{\tau }{\int }}{S}_i^{\pm }(t)d{S}_{\ell }(t)+\frac{S_i(\tau )S_{\ell }(\tau )}{2}\end{array}$$

The last term on the right-hand side corrects the reduced amount of the pairwise effect due to the truncated integral. Note that $w_{ii}=1/2$ . Accordingly, the estimator is also modified as

$$\begin{array}{r}{\hat{w}}_{i\ell }^{\text{CG}}=-\underset{0}{\overset{\tau }{\int }}{\hat{S}}_i^{\text{CG}\pm }(t)d{\hat{S}}_{\ell }^{\text{CG}}(t)+\frac{{\hat{S}}_i^{\text{CG}}(\tau ){\hat{S}}_{\ell }^{\text{CG}}(\tau )}{2}\end{array}$$

3.5 Sensitivity analysis

The proposed method needs to specify a copula (including its parameter ${\theta }_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 ${\theta }_i$ .^{30,31,47–49} If treatment effects are always significant irrespective of the values of ${\theta }_i$ , this result confirms the presence of treatment effects under dependent censoring. On the other hand, if the significance depends on the values of ${\theta }_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. ${\theta }_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{\mathit{p}}}^{\text{CG}}$ , their SE, and 95% CI (Section 3.2). It also computes the F-statistic and its critical value for testing $H_0:\mathit{C}\mathit{p}=\mathbf{\text{0}}$ under factorial designs (Section 3.3). The contrast matrix $\mathit{C}$ and the follow-up end $\tau$ 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{\theta }_1}=\cdots =C_{d{\theta }_d}$ and ${\theta }_1=\cdots ={\theta }_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{\mathit{p}}}^{\text{CG}}$ and the test statistic $F_N^{\text{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_{ij},U_{ij})$ , $i=1,2,\dots ,d,j=1,\dots ,n_i$ from the Clayton copula model with ${\theta }_i=2$ (Kendall's tau = 0.5) and exponential marginals with the rate parameters ${\lambda }_i$ for $T_{ij}$ and ${\mu }_i$ for $U_{ij}$ . We redefined $U_{ij}$ by $U_{ij}=\text{min}(U_{ij},\tau )$ , where $\tau =1$ is the follow-up end. We then obtained the censored data $\{(X_{ij},{\delta }_{ij});i=1,\dots ,d,j=1,\dots ,n_i\}$ , where $X_{ij}=\text{min}(T_{ij},U_{ij})$ and ${\delta }_{ij}=1\{T_{ij}\le U_{ij}\}$ .

Based on the data, we calculated the estimator ${\hat{p}}_i^{\text{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^{\text{CG}}$ , for testing

$$\begin{array}{r}{H}_0:p_1=p_2=p_3\text{versus}{H}_1:p_1\ne p_2,\text{or}{p}_1\ne p_3,\text{or}{p}_2\ne p_3\end{array}$$

by setting the contrast matrix

$$\begin{array}{r}\mathit{C}={\mathit{I}}_3-\frac{{\mathbf{\text{1}}}_3{\mathbf{\text{1}}}_3^{\mathrm{\prime }}}{3}=\frac{1}{3}\left(\begin{array}{ccc}2& -1& -1\\ -1& 2& -1\\ -1& -1& 2\end{array}\right)\end{array}$$

For the two-way layout ( $2\times 3=6$ levels), the true treatment effects $\mathit{p}$ are

$$\begin{array}{r}\mathit{p}=(p_1,p_2,p_3,p_4,p_5,p_6{)}^{\mathrm{\prime }}={\left(\underset{\text{Firstlevel in FactorA}}{\underbrace{p_{11},p_{12},p_{13}}},\underset{\text{Secondlevel in Factor A}}{\underbrace{p_{21},p_{22},p_{23}}}\right)}^{\mathrm{\prime }}\end{array}$$

To test the null effect of Factor A, we calculated the F-statistic, $F_N^{\text{CG}}$ , for no main effect in Factor A,

$$\begin{array}{r}{H}_0:p_{11}+p_{12}+p_{13}=p_{21}+p_{22}+p_{23}\text{versus}{H}_1:p_{11}+p_{12}+p_{13}\ne p_{21}+p_{22}+p_{23}\end{array}$$

by setting

$$\begin{array}{r}\mathit{C}={\mathit{P}}_2\otimes \frac{{\mathbf{\text{1}}}_3^{\mathrm{\prime }}}{3}=\frac{1}{6}\left(\begin{array}{cc}\begin{array}{ccc}1& 1& 1\\ -1& -1& -1\end{array}& \begin{array}{ccc}-1& -1& -1\\ 1& 1& 1\end{array}\end{array}\right)\end{array}$$

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

The rate parameters $({\lambda }_i,{\mu }_i)$ in the data-generating process were set to yield a realistic amount of treatment effects and censoring percentage $P(U_{ij}<T_{ij})\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_{ij}$	Rate parameters for $U_{ij}$	True treatment effects
One-way	$({\lambda }_1,{\lambda }_2,{\lambda }_3)$	$({\mu }_1,{\mu }_2,{\mu }_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	$({\lambda }_{11},{\lambda }_{12},{\lambda }_{13},{\lambda }_{21},{\lambda }_{22},{\lambda }_{23})$	$({\mu }_{11},{\mu }_{12},{\mu }_{13},{\mu }_{21},{\mu }_{22},{\mu }_{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 $\mathit{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 ${\theta }_i=2$ is known. Next, we turn our attention to the case, where ${\theta }_i=2$ is possibly misspecified as ${\theta }_i=\text{0, 1, 2, 3, and 4}$ .

When the true value ${\theta }_i=2$ is known, the means of ${\hat{p}}_i^{\text{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^{\text{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 ${\theta }_i=2$ is known. The type I error rates are well-controlled for the nominal level of $\alpha \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 $\alpha$ .

				Simulation method a			Analytical method b
	Design	Scenario	$n_i$	$\alpha =0.10$	$\alpha =0.05$	$\alpha =0.01$	$\alpha =0.10$	$\alpha =0.05$	$\alpha =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

^a

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

^b

The critical value is $c_{N,\alpha }={\chi }_{\hat{f},1-\alpha }^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^{\text{CG}}$ under the misspecification of ${\theta }_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 ${\theta }_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 ${\theta }_i$ deviates from the true value ${\theta }_i=2$ . Scenario 3 also yields some biases caused by misspecified values of ${\theta }_i$ . However, even when ${\theta }_i$ is misspecified and the estimator ${\hat{p}}_i^{\text{CG}}$ is biased, the coverage probabilities remain close to 0.95 in many cases. However, for the misspecified value of ${\theta }_i=0$ , the obvious under-coverage is found.

Figure 1.

Simulation results for ${\hat{p}}_i^{\text{CG}},i=1,2,3$ , when ${\theta }_i=2$ (true value) is misspecified as ${\theta }_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,\text{and}{p}_3=0.455$ .

Therefore, if the independent censoring assumption is mistakenly imposed, the estimator ${\hat{p}}_i^{\text{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^{\text{CG}},i=1,2,3$ , when ${\theta }_i=2$ (true value) is misspecified as ${\theta }_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, ${\theta }_i=0$ , ${\theta }_i=2$ , ${\theta }_i=4$ , and ${\theta }_i=8$ . Another obvious pattern is that estimated survival probabilities reduce remarkably when the copula parameter moves from ${\theta }_i=0$ to ${\theta }_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 ${\theta }_i=0$ , ${\theta }_i=2$ , ${\theta }_i=4$ , and ${\theta }_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)\mathrm{\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 ${\theta }_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 ${\theta }_i=0$ , ${\theta }_i=2$ , ${\theta }_i=4$ , and ${\theta }_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 ${\theta }_i$ . Hence, even if the independent censoring assumption is violated, the hypothesis $H_0^S:S_1(t)=S_2(t)=S_3(t)\mathrm{\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^{\text{CG}}$	$c_{N,\alpha }$ : simulation	$c_{N,\alpha }$ : analytical	P-value
$\theta =0$ (Kendall's tau = 0.00)	133.979	3.073	3.044	<0.0001
$\theta =2$ (Kendall's tau = 0.50)	15.453	3.763	3.515	0.0009
$\theta =4$ (Kendall's tau = 0.67)	22.764	3.602	3.405	0.0001
$\theta =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 ${\theta }_i$ . For instance, if ${\theta }_i=8$ is assumed, the adjuvant treatment yields a survival advantage $({\hat{p}}_1^{\text{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 ${\theta }_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
$\theta =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)
$\theta =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)
$\theta =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)
$\theta =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,\tau ]$ . 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