# Two-Stage Pseudo Maximum Likelihood Estimation of Semiparametric Copula-
based Regression Models for Semi-Competing Risks Data

Sakie J. Arachchige
Department of Mathematics and Statistics,
Mississippi State University, MS 39762, USA  Xinyuan Chen
Department of Mathematics and Statistics,
Mississippi State University, MS 39762, USA  Qian M. Zhou∗
Department of Mathematics and Statistics,
Mississippi State University, MS 39762, USA
$*$email<EMAIL_ADDRESS>
###### Abstract

We propose a two-stage estimation procedure for a copula-based model with
semi-competing risks data, where the non-terminal event is subject to
dependent censoring by the terminal event, and both events are subject to
independent censoring. Under a copula-based model, the marginal survival
functions of individual event times are specified by semiparametric
transformation models, and the dependence between the bivariate event times is
specified by a parametric copula function. For the estimation procedure, in
the first stage, the parameters associated with the marginal of the terminal
event are estimated only using the corresponding observed outcomes, and in the
second stage, the marginal parameters for the non-terminal event time and the
copula parameter are estimated via maximizing a pseudo-likelihood function
based on the joint distribution of the bivariate event times. We derived the
asymptotic properties of the proposed estimator and provided an analytic
variance estimator for inference. Through simulation studies, we showed that
our approach leads to consistent estimates with less computational cost and
more robustness compared to the one-stage procedure developed in Chen, (2012),
where all parameters were estimated simultaneously. In addition, our approach
demonstrates more desirable finite-sample performances over another existing
two-stage estimation method proposed in Zhu et al., (2021).

Key words: Copulas; Pseudo maximum likelihood estimator; Semi-competing risks;
Semiparametric regression; Stage-wise estimation.

##  1 Introduction

Medical studies oftentimes involve situations where multiple time-to-event
outcomes exist simultaneously. These events can be divided into two
categories: _non-terminal_ , e.g., cancer relapse, and _terminal_ , e.g.,
death. The occurrence of the terminal event precludes that of the non-
terminal, but not vice versa; the resulting data are commonly referred to as
the _semi-competing risks_ (Fine et al., 2001). Let $T$ and $D$ denote times
to a non-terminal and terminal event, respectively, and they are assumed to be
subject to independent administrative censoring. The distribution of $T$ and
the dependence between $T$ and $D$ are often of interest.

Copula-based models have been used in analyzing semi-competing risks data
(Fine et al., 2001, Wang, 2003, Lakhal et al., 2008), where the joint survival
function of $(T,D)$, $S(t,d)\equiv\mathbb{P}(T>t,D>d)$, is directly specified
as $S(t,d)=\mathbbm{C}(S_{T}(t),S_{D}(d))$. Here,
$\mathbbm{C}(\cdot,\cdot):[0,1]^{2}\rightarrow[0,1]$ is a bivariate copula
function (Sklar, 1959), and $S_{T}(t)\equiv\mathbb{P}(T>t)$ and
$S_{D}(d)\equiv\mathbb{P}(D>d)$ are marginal survival functions of $T$ and
$D$, respectively. This copula-based approach enjoys the flexibility from the
separate modeling of (i) the marginals, and (ii) the dependence structure
between $T$ and $D$. Shih and Louis, (1995) proposed a two-stage estimation
procedure for multivariate survival outcome under independent censoring,
where, in the first stage, the marginals were estimated separately for each
event without accounting for the between-event dependencies, and, in the
second stage, the parameter associated with the copula was estimated via
maximizing the log-likelihood with the plug-in estimates of the marginals. The
resulting estimator of the copula parameter is referred to as the pseudo
maximum likelihood estimator (PMLE).

An accurate estimation of the copula parameter relies on consistent estimation
of the marginals. The validity of the above stage-wise estimation procedure
requires that individual event times are only subject to independent
censoring. Thus, such an approach has been mostly used on clustered survival
data. However, this assumption does not hold for semi-competing risks data
because $T$ is dependently censored by $D$. In other words, standard tools
such as the Kaplan-Meier (KM) estimator are not consistent for $S_{T}(\cdot)$,
which poses challenges in estimating the copula parameter via the above two-
stage procedure. As a result, earlier estimators were obtained from estimating
equations based on a concordance indicator (Fine et al., 2001), the Doob–Meyer
(DM) decomposition of a counting process (Wang, 2003), or conditional
Kendall’s $\tau$ (Lakhal et al., 2008). Later, these procedures were extended
to incorporate a discrete covariate for the marginal distribution via a
regression model and for the between-event dependence parameter (Ghosh, 2006,
Hsieh et al., 2008). For estimating the marginal of $T$, they predominantly
considered a _plug-in_ estimator, first proposed in Fine et al., (2001),
through the relationship of $S_{T}(\cdot)$ with $S_{D}(\cdot)$ and the
survival function $S_{T^{*}}(\cdot)$ of $T^{*}\equiv T\wedge D$, the time to
the first event (non-terminal or terminal). However, this method is limited to
the Archimedean copula family, which provides an explicit expression for such
a relationship. Different from the above methods, Peng and Fine, (2007)
considered a class of time-varying coefficient semiparametric regression
models for the marginals and copulas with time-varying parameters. They
estimated the copula parameter and the marginal of $T$ simultaneously using
non-linear estimating equations.

Chen, (2012) considered the maximum likelihood estimation (MLE) approach,
i.e., simultaneous estimation of all parameters that are associated with the
marginals and the copula via maximizing the log-likelihood function. Zhou et
al., (2016) proposed a PMLE for the copula parameter based on the log-
likelihood function expressed in terms of $S_{D}(\cdot)$ and
$S_{T^{*}}(\cdot)$ without covariates. Zhu et al., (2021) extended the
procedure in Zhou et al., (2016) to incorporate covariates into the marginals
specified by time-varying coefficient semiparametric regression models, and
proposed a two-stage estimation procedure: in the first stage, the PMLE of the
copula parameter was obtained following the idea in Zhou et al., (2016), and,
in the second stage, the regression coefficients associated with
$S_{T}(\cdot)$ were estimated from non-linear estimating equations developed
in Peng and Fine, (2007).

In this article, we propose a new estimation procedure for copula-based models
to analyze semi-competing risks data, with the primary goals of investigating
(i) the association between $T$ and $D$, and (ii) the direct covariate effects
on their marginals. Along the line of Chen, (2012), we consider a class of
models where a parametric copula is employed to model the dependence between
$T$ and $D$, and the marginals are specified by semi-parametric transformation
models. We estimate the marginals and copula parameters via a two-stage PMLE
approach. Specifically, in the first stage, we obtain a consistent estimate of
the marginal of $D$, since it is only subject to independent censoring, and,
in the second stage, we estimate the marginal of $T$ and the copula parameter
simultaneously by maximizing the pseudo-likelihood function with the plug-in
estimates of $S_{D}(\cdot)$.

Our two-stage estimation approach adapts from the one-stage MLE method from
Chen, (2012), which is expected to be more efficient but likely
computationally challenging in practice. In contrast, our PMLE approach
reduces the dimension of the parameter space in each stage while yielding
consistent estimates. Despite not utilizing the association between $T$ and
$D$ in the first stage, which might result in efficiency loss, we observed
through simulation studies that, compared to the MLE, the PMLE approach gives
more desirable finite-sample performances including smaller biases for the
marginal parameters of $T$ and copula parameters. In addition, our approach
enjoys robustness against copula misspecification for estimating the
marginals. When the copula function is misspecified, our estimation of
$S_{D}(\cdot)$ in the first stage is guaranteed to be consistent. Furthermore,
although the PMLE of $S_{T}(\cdot)$ is not consistent under copula
misspecification, with a consistent estimate of $S_{D}(\cdot)$, it leads to a
smaller bias compared to the one-stage MLE, which is additionally affected by
the biased estimates of $S_{D}(\cdot)$.

Compared to the stage-wise estimation procedure developed by Zhu et al.,
(2021), our method exhibits three advantages. First, our method does not
involve the first-event time $T^{*}$, and consequently, besides Archimedean
copulas, it can also be applied to the elliptical copula family, e.g.,
Gaussian or $t$ copulas. Second, our estimators for the marginal of $T$ and
the copula parameter are inherently more efficient, since, assuming the
marginal of $D$ is known, our estimator is equivalent to the MLE, but the
estimator in Zhu et al., (2021) is still the PMLE based on plug-in estimates.
Third, we provide variance estimators based on analytic expressions for
inference. Compared to the non-parametric bootstrap employed in Zhu et al.,
(2021) and Zhou et al., (2016), our inference procedure is computationally
more efficient. To implement our method, an R package PMLE4SCR is available in
the GitHub repository https://github.com/michellezhou2009/PMLE4SCR.

The remainder of the article is organized as follows. In Section 2, we propose
the two-stage PMLE estimation procedure and discuss its asymptotic properties.
Section 3 presents simulation studies evaluating the finite-sample performance
of the PMLE and the analysis of Bone Marrow Transplant (BMT) data as an
illustration. We provide concluding remarks in Section 4.

##  2 Method

###  2.1 Copula-based model specification

Let $A$ denote the time to independent administrative censoring. The observed
event times are denoted by $C\equiv D\wedge A$ and $X\equiv T\wedge C$, with
respective censoring status defined as $\delta_{D}=\mathbbm{1}(D\leq A)$ and
$\delta_{T}=\mathbbm{1}(T\leq C)$. In addition, let ${\bf Z}\in\mathbb{R}^{p}$
denote a vector of baseline covariates, and given ${\bf Z}$, $A$ is
independent of both $T$ and $D$, but $T$ and $D$ are associated. The observed
data from $n$ independent subjects is denoted by
$\mathfrak{S}=\\{(X_{i},C_{i},\delta_{T,i},\delta_{D,i},{\bf
Z}_{i}),i=1,\ldots,n\\}$. We have the following commonly adopted assumptions
regarding the study setting and covariates:

  * (A1)

There exists a maximum follow-up time $\xi<\infty$ in the study such that
$\mathbb{P}(X>\xi|{\bf Z})>0$,
$\mathbb{P}(\delta_{T}=0,\delta_{D}=0,X=\xi|{\bf Z})>0$,
$\mathbb{P}(C>\xi|{\bf Z})>0$, and
$\mathbb{P}(\delta_{T}=1,\delta_{D}=0,C=\xi|{\bf Z})>0$, with probability 1;

  * (A2)

The covariates ${\bf Z}$ is bounded with probability 1.

Given covariates ${\bf Z}$, the copula-based approach models the joint
survival function of $(T,D)$, $\mathbb{P}(T>t,D>d|{\bf Z})$, as a function of
their respective marginals $S_{T}(t|{\bf Z})=\mathbb{P}(T>t|{\bf Z})$ and
$S_{D}(d|{\bf Z})=\mathbb{P}(D>d|{\bf Z})$, i.e.,

$\mathbb{P}(T>t,D>d|{\bf Z})=\mathbbm{C}\left(S_{T}(t|{\bf Z}),S_{D}(d|{\bf
Z});\alpha\right),$ (1)

where $\mathbbm{C}(\cdot,\cdot;\alpha):[0,1]^{2}\rightarrow[0,1]$ is a copula
function parametrically specified with parameter $\alpha$. In the remainder of
this article, we define
$\dot{\mathbbm{C}}_{u_{j}}\equiv\partial\mathbbm{C}(u_{T},u_{D};\alpha)/\partial
u_{j}$, $j=T,D$, as the first-order partial derivatives of the copula
function, and
$\ddot{\mathbbm{C}}_{u_{T}u_{D}}\equiv\partial^{2}\mathbbm{C}(u_{T},u_{D};\alpha)/\partial
u_{T}u_{D}$ as the second-order partial derivative.

###### Remark 1

The Archimedean and Elliptical copula families are often selected for
specifying $\mathbbm{C}(\cdot,\cdot;\alpha)$. Archimedean copulas are defined
via a generator function $\zeta(u;\alpha)$ for $u\in[0,1]$, which is
continuous, convex, and strictly decreasing with $\zeta(1;\alpha)=0$. A
bivariate Archimedean copula can be constructed as
$\mathbbm{C}(u_{1},u_{2};\alpha)=\zeta^{-1}\\{\zeta(u_{1};\alpha)+\zeta(u_{2};\alpha);\alpha\\}$.
We present some commonly used generator functions and corresponding bivariate
copulas below and refer readers to Nelsen, (2006) for a more detailed
introduction.

  1. (i)

Clayton copula, $\zeta(u;\alpha)=(u^{-\alpha}-1)/\alpha$ for
$\alpha\in(-1,\infty)\backslash\\{0\\}$, and

$\displaystyle\mathbbm{C}(u_{1},u_{2};\alpha)=\left\\{\max\left(u_{1}^{-\alpha}+u_{2}^{-\alpha}-1,0\right)\right\\}^{-1/\alpha};$
(2)

  2. (ii)

Frank copula, $\zeta(u;\alpha)=-\log\\{(e^{-\alpha u}-1)/(e^{-\alpha}-1)\\}$
for $\alpha\in\mathbb{R}\backslash\\{0\\}$, and

$\displaystyle\mathbbm{C}(u_{1},u_{2};\alpha)=-\alpha^{-1}\log\left\\{1+\frac{(e^{-\alpha
u_{1}}-1)(e^{-\alpha u_{2}}-1)}{e^{-\alpha}-1}\right\\};$ (3)

  3. (iii)

Gumbel copula, $\zeta(u;\alpha)=\\{-\log(u)\\}^{\alpha}$ for
$\alpha\in[1,\infty)$, and

$\displaystyle\mathbbm{C}(u_{1},u_{2};\alpha)=\exp\left(-\left[\\{-\log(u_{1})\\}^{\alpha}+\\{-\log(u_{2})\\}^{\alpha}\right]^{1/\alpha}\right).$
(4)

Elliptical copulas join marginals via an elliptical distribution, e.g., a
Gaussian distribution. A bivariate Gaussian copula with parameter $\alpha$ is
defined as

$\displaystyle\mathbbm{C}(u_{1},u_{2};\alpha)=\Phi_{2}\\{\Phi^{-1}(u_{1}),\Phi^{-1}(u_{2});\alpha\\},$
(5)

where $\Phi(\cdot)$ is the univariate cumulative distribution function of the
standard Gaussian distribution, and $\Phi_{2}(\cdot,\cdot;\alpha)$ is the
standard bivariate Gaussian cumulative distribution function with $\alpha$
being the correlation parameter.

We specify the marginal survival functions of $T$ and $D$ via semi-parametric
transformation models. Specifically, for $j=T,D$,

$S_{j}(t|{\bf
Z})=\exp\left[-G_{j}\left\\{R_{j}(t)e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}}\right\\}\right],$ (6)

where $\boldsymbol{\beta}_{j}$ is a $p$-dimensional vector of regression
coefficients, the unspecified baseline function $R_{j}(\cdot)$ is assumed to
be increasing, and $G_{j}(\cdot)$ is a specified non-negative, strictly
increasing, and continuously differentiable function. For example, if
$G_{j}(t)=t$, (6) becomes a Cox proportional hazards (PH) model; if
$G_{j}(t)=\log(1+t)$, (6) becomes a proportional odds model. The marginal
survival function in (6) can be easily extended to accommodate time-varying
covariates ${\bf Z}(t)$, where $S_{j}(t|{\bf
Z})=\exp[-G_{j}\\{R_{j}(t)e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}(t)}\\}]$.

###  2.2 The two-stage estimation procedure

Let $\boldsymbol{\theta}_{D}=\\{\boldsymbol{\beta}_{D},R_{D}(\cdot)\\}$ and
$\boldsymbol{\theta}_{T}=\\{\boldsymbol{\beta}_{T},R_{T}(\cdot)\\}$ denote the
parameters associated with the marginals of $D$ and $T$, respectively. We
propose a two-stage estimation procedure for
$\boldsymbol{\theta}=\\{\boldsymbol{\theta}_{D},\boldsymbol{\theta}_{T},\alpha\\}$,
where, in the first stage, $\boldsymbol{\theta}_{D}$ is estimated separately
via an MLE approach with observed data
$\mathfrak{S}_{D}=\\{(C_{i},\delta_{D,i},{\bf Z}_{i}),i=1,\ldots,n\\}$, and,
in the second stage,
$\boldsymbol{\theta}_{1}\equiv\\{\boldsymbol{\theta}_{T},\alpha\\}$ is
estimated via maximizing a pseudo-log-likelihood depending on the estimate
$\widehat{\boldsymbol{\theta}}_{D}$ obtained from the first stage. We
hereafter refer to the proposed estimation procedure as the two-stage PMLE for
brevity.

####  2.2.1 The first stage estimation of $\boldsymbol{\theta}_{D}$

In the semi-competing risks setting, $\boldsymbol{\theta}_{D}$ can be
consistently estimated by an MLE with observed data $\mathfrak{S}_{D}$, since
$D$ is only subject to independent censoring by $A$ (Zeng and Lin, 2006). The
baseline function $R_{D}(\cdot)$ is regarded as a non-decreasing step function
that jumps only at the observed event times. Let
$0<d_{1}<d_{2}<\cdots<d_{\kappa_{D}}\leq\xi$ denote the _ordered_ observed
death times, and let $\mathrm{d}R_{D,l}>0$ denote the jump size of
$R_{D}(\cdot)$ at time $d_{l}$, for $l=1,\ldots,\kappa_{D}$. Thus,
$\boldsymbol{\theta}_{D}$ can be re-expressed as
$\boldsymbol{\theta}_{D}=(\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}},\mathrm{d}R_{D,1},\ldots,\mathrm{d}R_{D,\kappa_{D}})^{{\sf\scriptscriptstyle{T}}}$.
The log-likelihood function of $\boldsymbol{\theta}_{D}$ given
$\mathfrak{S}_{D}$ is
$\ell_{D}(\boldsymbol{\theta}_{D})=n^{-1}\sum_{i=1}^{n}\ell_{D,i}(\boldsymbol{\theta}_{D})$,
where

$\ell_{D,i}(\boldsymbol{\theta}_{D})=\int_{0}^{\xi}\left[\log\mathrm{d}R_{D}(t)+\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}+\log\dot{G}_{D}\left\\{R_{D}(t)e^{\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\right\\}\right]\mathrm{d}N_{D,i}(t)-G_{D}\left\\{R_{D}(C_{i})e^{\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\right\\},$ (7)

with $N_{D,i}(t)=\mathbbm{1}(C_{i}\leq t)\delta_{D,i}$. The MLE
$\widehat{\boldsymbol{\theta}}_{D}=(\widehat{\boldsymbol{\beta}}_{D}^{{\sf\scriptscriptstyle{T}}},\widehat{\mathrm{d}R}_{D,1},\ldots,\widehat{\mathrm{d}R}_{D,\kappa_{D}})^{{\sf\scriptscriptstyle{T}}}$
is the solution to score equations
$\Psi_{D}(\boldsymbol{\theta}_{D})=n^{-1}\sum_{i=1}^{n}\Psi_{D,i}(\boldsymbol{\theta}_{D})=0$
with
$\Psi_{D,i}(\boldsymbol{\theta}_{D})=\partial\ell_{D,i}(\boldsymbol{\theta}_{D})/\partial\boldsymbol{\theta}_{D}$.

####  2.2.2 The second stage estimation of $\boldsymbol{\theta}_{1}$

Different from $\boldsymbol{\theta}_{D}$, $\boldsymbol{\theta}_{T}$ cannot be
consistently estimated using only the observed data
$\\{(X_{i},\delta_{T,i},{\bf Z}_{i}),\allowbreak i=1,\ldots,n\\}$, since $T$
is subject to both dependent and independent censoring by $D$ and $A$,
respectively. A consistent estimator must be obtained from a full likelihood
function based on the joint survival function of $(T,D)$ with $\mathfrak{S}$
(Fine et al., 2001). Similar to $R_{D}(\cdot)$, $R_{T}(\cdot)$ is also
regarded as a non-decreasing step function that jumps only at the observed
non-terminal event times $0<t_{1}<t_{2}<\cdots<t_{\kappa_{T}}\leq\xi$, where
$\mathrm{d}R_{T,l}$ denotes the jump size of $R_{T}(\cdot)$ at time $t_{l}$,
for $l=1,\ldots,\kappa_{T}$. The parameter $\boldsymbol{\theta}_{T}$ can be
re-expressed as
$\boldsymbol{\theta}_{T}=(\boldsymbol{\beta}_{T}^{{\sf\scriptscriptstyle{T}}},\mathrm{d}R_{T,1},\ldots,\mathrm{d}R_{T,\kappa_{T}})^{{\sf\scriptscriptstyle{T}}}$.
The full log-likelihood function given $\mathfrak{S}$ is
$\ell(\boldsymbol{\theta})=n^{-1}\sum_{i=1}^{n}\ell_{i}(\boldsymbol{\theta})$,
where

$\displaystyle\ell_{i}(\boldsymbol{\theta})$
$\displaystyle=\delta_{T,i}\delta_{D,i}\log\ddot{\mathbbm{C}}_{u_{T}u_{D}}+\delta_{T,i}(1-\delta_{D,i})\log\dot{\mathbbm{C}}_{u_{T}}$
(8)
$\displaystyle~{}+(1-\delta_{T,i})\delta_{D,i}\log\dot{\mathbbm{C}}_{u_{D}}+(1-\delta_{T,i})(1-\delta_{D,i})\log\mathbbm{C}$
$\displaystyle~{}+\delta_{T,i}\left[-G_{T}\left\\{R_{T}(X_{i})e^{\boldsymbol{\beta}_{T}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\right\\}+\log\dot{G}_{T}\left\\{R_{T}(X_{i})e^{\boldsymbol{\beta}_{T}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\right\\}+\log\mathrm{d}R_{T}(X_{i})+\boldsymbol{\beta}_{T}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}\right]$
$\displaystyle~{}+\delta_{D,i}\left[-G_{D}\left\\{R_{D}(C_{i})e^{\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\right\\}+\log\dot{G}_{D}\left\\{R_{D}(C_{i})e^{\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\right\\}+\log\mathrm{d}R_{D}(C_{i})+\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}\right],$

with
$U_{T,i}=\exp\left[-G_{T}\left\\{R_{T}(X_{i})e^{\boldsymbol{\beta}_{T}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\right\\}\right]$ and
$U_{D,i}=\exp\left[-G_{D}\left\\{R_{D}(C_{i})e^{\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\right\\}\right]$. The PMLE of $\boldsymbol{\theta}_{1}$ is obtained by
maximizing the following pseudo-log-likelihood function:

$\displaystyle\widehat{\ell}(\boldsymbol{\theta}_{1})\equiv\ell(\boldsymbol{\theta}_{1},\widehat{\boldsymbol{\theta}}_{D})=n^{-1}\sum_{i=1}^{n}\ell_{i}(\boldsymbol{\theta}_{1},\widehat{\boldsymbol{\theta}}_{D}),$

where $\ell_{i}$ is given in (8). The PMLE
$\widehat{\boldsymbol{\theta}}_{1}=(\widehat{\alpha},\widehat{\boldsymbol{\beta}}_{T}^{{\sf\scriptscriptstyle{T}}},\widehat{\mathrm{d}R}_{T,1},\ldots,\widehat{\mathrm{d}R}_{T,\kappa_{T}})^{{\sf\scriptscriptstyle{T}}}$
is the solution to pseudo score equations
$\widehat{\Psi}_{1}(\boldsymbol{\theta}_{1})=n^{-1}\sum_{i=1}^{n}\widehat{\Psi}_{1,i}(\boldsymbol{\theta}_{1})=0$
with
$\widehat{\Psi}_{1,i}(\boldsymbol{\theta}_{1})=\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1},\widehat{\boldsymbol{\theta}}_{D})=\partial\ell_{i}(\boldsymbol{\theta}_{1},\widehat{\boldsymbol{\theta}}_{D})/\partial\boldsymbol{\theta}_{1}$.
The score function
$\ell_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{D})=(\Psi_{\alpha,i},\Psi_{\boldsymbol{\beta}_{T},i}^{{\sf\scriptscriptstyle{T}}},\Psi_{\mathrm{d}R_{T1},i},\ldots,\Psi_{\mathrm{d}R_{T\kappa_{T}},i})^{{\sf\scriptscriptstyle{T}}}$,
with the expressions given in Appendix A, is a function of
$\boldsymbol{\theta}_{D}$ only through $U_{D,i}$.

###### Remark 2

The copula in (1) can be extended by allowing the copula parameter to depend
on some baseline covariates. For example, we can assume
$\alpha=\phi(\boldsymbol{\gamma}^{{\sf\scriptscriptstyle{T}}}{\bf W})$, where
${\bf W}$ is a subset of the covariates ${\bf Z}$, and $\phi(\cdot)$ is a
specified link function. The copula parameter thus becomes
$\boldsymbol{\gamma}$, and its corresponding score function is modified as
$\Psi_{\boldsymbol{\gamma},i}=\Psi_{\alpha,i}\dot{\phi}(\boldsymbol{\gamma}^{{\sf\scriptscriptstyle{T}}}{\bf
W}_{i}){\bf W}_{i}$, where $\dot{f}(\cdot)$ denotes the first-order derivative
of a univariate function $f(\cdot)$.

###### Remark 3

The log-likelihood function in (8) is equivalent to the log-likelihood in
Chen, (2012) derived based on three observed counting processes: the first
event being non-terminal, the first event being terminal, and the second event
occurring subsequently to the uncensored first event. Therefore, the one-stage
MLE obtained via maximizing (8) shares the same asymptotic properties given in
Theorems 1 and 2 of Chen, (2012).

###  2.3 Asymptotic properties of the two-stage PMLE

We present the asymptotic properties of the two-stage PMLE. Define the
following quantities:

$\displaystyle\Upsilon_{1}(g_{T},g_{D})=-\log\mathbbm{C}\left(e^{-G_{T}(g_{T})},e^{-G_{D}(g_{D})};\alpha\right)~{}~{}\text{and}~{}~{}\eta_{1j}(g_{T},g_{D})=\frac{\partial}{\partial
g_{j}}\Upsilon_{1}(g_{T},g_{D}),$

for $j=T,D$; also

$\displaystyle\Upsilon_{2}(g_{T},g_{D})=-\log\dot{\mathbbm{C}}_{u_{T}}\left(e^{-G_{T}(g_{T})},e^{-G_{D}(g_{D})};\alpha\right)~{}~{}\text{and}~{}~{}\eta_{2}(g_{T},g_{D})=\frac{\partial}{\partial
g_{D}}\Upsilon_{2}(g_{T},g_{D}).$

Let
$\boldsymbol{\theta}^{0}=\\{\alpha^{0},\boldsymbol{\beta}_{T}^{0},R_{T}^{0}(\cdot),\allowbreak\boldsymbol{\beta}_{D}^{0},R_{D}^{0}(\cdot)\\}$
be the true values of $\boldsymbol{\theta}$, and $\boldsymbol{\Theta}$ denote
the total parameter space. The following conditions are required to establish
the asymptotic properties of
$\widehat{\boldsymbol{\theta}}=(\widehat{\boldsymbol{\theta}}_{D}^{{\sf\scriptscriptstyle{T}}},\widehat{\boldsymbol{\theta}}_{1}^{{\sf\scriptscriptstyle{T}}})^{{\sf\scriptscriptstyle{T}}}$.

  * (C1)

The true baseline functions $R_{T}^{0}(\cdot)$ and $R_{D}^{0}(\cdot)$ are
strictly increasing and continuously differentiable, and the true values
$\alpha^{0}$, $\boldsymbol{\beta}_{T}^{0}$ and $\boldsymbol{\beta}_{D}^{0}$,
fall in the interior of a compact parameter space.

  * (C2)

For any sequence $0<h_{1}<\cdots<h_{r}\leq g$,

$\prod_{s=1}^{r}\\{(1+h_{s})\dot{G}_{D}(h_{s})\\}\exp\\{-G_{D}(g)\\}\leq\mu_{0}^{r}(1+g)^{-\nu_{0}},$

where $\mu_{0}$ and $\nu_{0}$ are positive constants.

  * (C3)

The function $G_{j}(\cdot)$ satisfies the conditions: $G_{j}(0)=0$,
$G_{j}(\infty)=\infty$, and $\dot{G}_{j}(0)>0$ for $j=T,D$. The copula
function $\mathbbm{C}(\cdot,\cdot;\alpha)$ satisfies the conditions:
$\mathbbm{C}(0,0;\alpha)=0$ and $\eta_{1j}(0,0)>0$ for $j=T,D$. Also, for
positive constants $c_{1}$, $c_{2}$, and $c_{3}$,

$\displaystyle\lim_{g_{1},g_{2}\rightarrow\infty}\frac{\sum_{j=1}^{2}\log\\{g_{j}\times\sup_{s_{1}\leq
g_{1},s_{2}\leq
g_{2}}\eta_{1j}(s_{1},s_{2})\\}}{\Upsilon_{1}(c_{1}g_{1},c_{2}g_{2})}=0,$
$\displaystyle\lim_{g_{1},g_{2}\rightarrow\infty}\frac{\log\\{g_{2}\times\sup_{s_{1}\leq
g_{1},s_{2}\leq g_{2}}\eta_{2}(s_{1},s_{2})\\}}{\inf_{w_{1}\geq 0,w_{2}\geq
0}\\{\Upsilon_{2}(c_{1}w_{1},c_{3}g_{2})-\Upsilon_{2}(c_{1}w_{1},c_{2}w_{2})\\}}=0.$

  * (C4)

$\sum_{i=1}^{n}\ell_{i}(\boldsymbol{\theta})$,
$\sum_{i=1}^{n}\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta})$,
$\sum_{i=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}_{1}\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta})$,
and
$\sum_{i=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},i}(\boldsymbol{\theta})$
are continuous and bounded on $[0,\xi]\times\boldsymbol{\Theta}$, where
$\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta})=\partial\ell_{i}(\boldsymbol{\theta})/\partial\boldsymbol{\theta}_{1}$,
$\ddot{\ell}_{\boldsymbol{\theta}_{1}\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta})=\partial^{2}\ell_{i}(\boldsymbol{\theta})/\partial\boldsymbol{\theta}_{1}\partial\boldsymbol{\theta}_{1}^{{\sf\scriptscriptstyle{T}}}$,
and
$\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},i}(\boldsymbol{\theta})=\partial^{2}\ell_{i}(\boldsymbol{\theta})/\partial\boldsymbol{\theta}_{1}\partial
u_{D}$.

  * (C5)

Information matrices
$\mathcal{I}_{D}(\boldsymbol{\theta}_{D})=\mathbb{E}\\{\mathcal{I}_{D,i}(\boldsymbol{\theta}_{D})\\}$
with
$\mathcal{I}_{D,i}(\boldsymbol{\theta}_{D})=-\partial^{2}\ell_{D,i}(\boldsymbol{\theta}_{D})/\partial\boldsymbol{\theta}_{D}\partial\boldsymbol{\theta}_{D}^{{\sf\scriptscriptstyle{T}}}$
and
$\mathcal{I}_{1}(\boldsymbol{\theta})=\mathbb{E}\\{\mathcal{I}_{1,i}(\boldsymbol{\theta})\\}$
with
$\mathcal{I}_{1,i}(\boldsymbol{\theta})=-\partial^{2}\ell_{i}(\boldsymbol{\theta})/\partial\boldsymbol{\theta}_{1}\partial\boldsymbol{\theta}_{1}^{{\sf\scriptscriptstyle{T}}}$
evaluated at the true value of $\boldsymbol{\theta}_{D}^{0}$ and
$\boldsymbol{\theta}^{0}$ are positive definite, and their eigenvalues are
bounded below and above by some finite positive constants. The explicit
expressions of $\mathcal{I}_{D}(\boldsymbol{\theta}_{D})$ and
$\mathcal{I}_{1}(\boldsymbol{\theta})$ are given in Appendix B.

  * (C6)

The variance-covariance matrix
$\mathrm{Var}\\{\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta})+n^{-1}\sum_{k=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},k}(\boldsymbol{\theta})\psi_{u_{D,k},i}(\boldsymbol{\theta}_{D})\\}$
evaluated at $\boldsymbol{\theta}^{0}$ is positive definite, and its
eigenvalues are bounded below and above by some finite positive constants,
where $\psi_{u_{D,k},i}(\boldsymbol{\theta}_{D})$ is given in Appendix C.

A1 - A2 and C1 - C2 are from the conditions listed in Zeng and Lin, (2006) for
establishing the consistency and asymptotic distribution of
$\widehat{\boldsymbol{\theta}}_{D}$. Specifically, with probability 1,
$\widehat{\boldsymbol{\beta}}_{D}$ converges to $\boldsymbol{\beta}_{D}^{0}$
and
$\widehat{R}_{D}(\cdot)\equiv\\{\widehat{\mathrm{d}R}_{D,l},l=1,\ldots,\kappa_{D}\\}$
converges to $R_{D}^{0}(\cdot)$ uniformly in the interval $[0,\xi]$. In
addition,
$\sqrt{n}(\widehat{\boldsymbol{\theta}}_{D}-\boldsymbol{\theta}_{D})$
converges weakly to a zero-mean Gaussian process. Let
$\mathcal{I}_{D}^{0}=\mathbb{E}\\{\mathcal{I}_{D,i}(\boldsymbol{\theta}_{D}^{0})\\}$.
For a Hamadard differentiable functional $F_{D}(\boldsymbol{\theta}_{D})$ of
$\boldsymbol{\theta}_{D}$ (Chen, 2012), the functional delta method suggests
that
$\sqrt{n}\\{F_{D}(\widehat{\boldsymbol{\theta}}_{D})-F_{D}(\boldsymbol{\theta}_{D})\\}$
converges weakly to a zero-mean Gaussian process with the variance-covariance
matrix
$\dot{F}_{D}(\boldsymbol{\theta}_{D})^{{\sf\scriptscriptstyle{T}}}(\mathcal{I}_{D}^{0})^{-1}\dot{F}_{D}(\boldsymbol{\theta}_{D})$,
where $\dot{F}_{D}(\boldsymbol{\theta}_{D})$ is the derivative of
$F_{D}(\boldsymbol{\theta}_{D})$ with respect to $\boldsymbol{\theta}_{D}$.
The matrix $(\mathcal{I}_{D}^{0})^{-1}$ is estimated by
$\mathcal{I}_{D,n}^{-1}(\widehat{\boldsymbol{\theta}}_{D})$ where
$\mathcal{I}_{D,n}(\widehat{\boldsymbol{\theta}}_{D})=n^{-1}\sum_{i=1}^{n}\mathcal{I}_{D,i}(\widehat{\boldsymbol{\theta}}_{D})$.
Alternatively, one can estimate $(\mathcal{I}_{D}^{0})^{-1}$ via a robust
estimator, given by $\widehat{\boldsymbol{\Sigma}}_{D,n}\equiv
n^{-1}\sum_{i=1}^{n}\psi_{\boldsymbol{\theta}_{D},i}(\widehat{\boldsymbol{\theta}}_{D})\psi_{\boldsymbol{\theta}_{D},i}(\widehat{\boldsymbol{\theta}}_{D})^{{\sf\scriptscriptstyle{T}}}$
with
$\psi_{\boldsymbol{\theta}_{D},i}(\widehat{\boldsymbol{\theta}}_{D})=\mathcal{I}_{D,n}^{-1}(\widehat{\boldsymbol{\theta}}_{D})\Psi_{D,i}(\widehat{\boldsymbol{\theta}}_{D})$,
where $\Psi_{D,i}(\boldsymbol{\theta}_{D})$ is the score function for the log-
likelihood function $\ell_{D,i}$ discussed in Section 2.2.1.

C3 regulates the tail behavior of the estimated baseline functions, which was
adopted in Chen, (2012). C4 and C5 are standard conditions guaranteeing the
existence of the score functions, information matrices, and certain components
for the variance-covariance matrix of $\widehat{\boldsymbol{\theta}}_{1}$. C6
assumes the variance-covariance matrix of $\widehat{\boldsymbol{\theta}}_{1}$
is well-defined. Additionally, since, for the Gumbel copula,
$\eta_{1j}(0,0)>0$ (in C3) does not hold, we adopt the same modification as
Chen, (2012) by adding a small positive value, e.g., $n^{-1}$, to the
cumulative hazard function,
$G_{j}\\{R_{j}(t)e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}}\\}$, at time 0 to make $\eta_{1j}(0,0)$ positive for $j=T,D$.

The consistency and asymptotic distribution of
$\widehat{\boldsymbol{\theta}}_{1}$ are established in the following Theorems
1 and 2, respectively.

###### Theorem 1

Assume A1 - A2 and C1 - C6 hold. With probability 1, $\widehat{\alpha}$
converges to $\alpha^{0}$, $\widehat{\boldsymbol{\beta}}_{T}$ converges to
$\boldsymbol{\beta}_{T}^{0}$, and
$\widehat{R}_{T}(\cdot)\equiv\\{\widehat{\mathrm{d}R}_{T,l},l=1,\ldots,\kappa_{T}\\}$
converges to $R_{T}^{0}(\cdot)$ uniformly in the interval $[0,\xi]$.

###### Theorem 2

Assume A1 - A2 and C1 - C6 hold, and then
$\sqrt{n}(\widehat{\boldsymbol{\theta}}_{1}-\boldsymbol{\theta}_{1})$
converges weakly to a zero-mean Gaussian process. Let
$\mathcal{I}_{1}^{0}=\mathbb{E}\\{\mathcal{I}_{1,i}(\boldsymbol{\theta}^{0})\\}$.
For a Hamadard differentiable functional $F_{1}(\boldsymbol{\theta}_{1})$ of
$\boldsymbol{\theta}_{1}$,
$\sqrt{n}\\{F_{1}(\widehat{\boldsymbol{\theta}}_{1})-F_{1}(\boldsymbol{\theta}_{1})\\}$
converges weakly to a zero-mean Gaussian process with the variance-covariance
matrix
$\dot{F}_{1}(\boldsymbol{\theta}_{1})^{{\sf\scriptscriptstyle{T}}}\boldsymbol{\Sigma}_{1}^{0}\dot{F}(\boldsymbol{\theta}_{1})$,
where $\dot{F}_{1}(\boldsymbol{\theta}_{1})$ is the derivative of
$F_{1}(\boldsymbol{\theta}_{1})$ with respect to $\boldsymbol{\theta}_{1}$,
and
$\boldsymbol{\Sigma}_{1}^{0}\equiv\mathrm{Var}\\{\psi_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}^{0})\\}$
with

$\displaystyle\psi_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}^{0})=\left(\mathcal{I}_{1}^{0}\right)^{-1}\left\\{\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}^{0})+n^{-1}\sum_{k=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},k}(\boldsymbol{\theta}^{0})\psi_{u_{D,k},i}(\boldsymbol{\theta}_{D}^{0})\right\\}.$

The expression of $\psi_{u_{D,k},i}(\boldsymbol{\theta}_{D}^{0})$ is given in
Appendix C.

An application of Theorem 2 is to consider linear functionals in the form of
$b_{\alpha}\widehat{\alpha}+{\bf
b}_{\beta_{T}}^{{\sf\scriptscriptstyle{T}}}\boldsymbol{\beta}_{T}+\int_{0}^{\xi}b_{R_{T}}(t)\mathrm{d}R_{T}(t)$
for $F_{1}(\boldsymbol{\theta}_{1})$, where $b_{\alpha}$ and ${\bf
b}_{\beta_{T}}$ are a known scalar and vector, respectively, and ${\bf
b}_{R_{T}}\equiv(b_{R_{T},1},\ldots,b_{R_{T},\kappa_{T}})^{{\sf\scriptscriptstyle{T}}}$
is the known vector of the values of $b_{R_{T}}(t)$ evaluated at
$0<t_{1}<t_{2}<\cdots<t_{\kappa_{T}}\leq\xi$. The corresponding variance-
covariance matrix is ${\bf
b}_{1}^{{\sf\scriptscriptstyle{T}}}\boldsymbol{\Sigma}_{1}^{0}{\bf b}_{1}$,
where ${\bf b}_{1}=(b_{\alpha},{\bf
b}_{\beta_{T}}^{{\sf\scriptscriptstyle{T}}},{\bf
b}_{R_{T}}^{{\sf\scriptscriptstyle{T}}})^{{\sf\scriptscriptstyle{T}}}$. The
variance-covariance matrix $\boldsymbol{\Sigma}_{1}^{0}$ can be consistently
estimated by

$\displaystyle\widehat{\boldsymbol{\Sigma}}_{1,n}=n^{-1}\sum_{i=1}^{n}\psi_{\boldsymbol{\theta}_{1},i}(\widehat{\boldsymbol{\theta}})\psi_{\boldsymbol{\theta}_{1},i}(\widehat{\boldsymbol{\theta}})^{{\sf\scriptscriptstyle{T}}},$

where

$\displaystyle\psi_{\boldsymbol{\theta}_{1},i}(\widehat{\boldsymbol{\theta}})=\mathcal{I}_{1,n}^{-1}(\widehat{\boldsymbol{\theta}})\left\\{\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\widehat{\boldsymbol{\theta}})+n^{-1}\sum_{k=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},k}(\widehat{\boldsymbol{\theta}})\psi_{u_{D,k},i}(\widehat{\boldsymbol{\theta}}_{D})\right\\}$

with
$\mathcal{I}_{1,n}(\widehat{\boldsymbol{\theta}})=n^{-1}\sum_{i=1}^{n}\mathcal{I}_{1,i}(\widehat{\boldsymbol{\theta}})$.
This estimator can be re-written as the form of a robust sandwich variance
estimator:
$\widehat{\boldsymbol{\Sigma}}_{1,n}=\mathcal{I}_{1,n}^{-1}(\widehat{\boldsymbol{\theta}})\mathcal{V}_{1,n}(\widehat{\boldsymbol{\theta}})\mathcal{I}_{1,n}^{-1}(\widehat{\boldsymbol{\theta}})$
with
$\mathcal{V}_{1,n}(\boldsymbol{\theta})=n^{-1}\sum_{i=1}^{n}[\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\widehat{\boldsymbol{\theta}})+n^{-1}\sum_{k=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},k}(\widehat{\boldsymbol{\theta}})\psi_{u_{D,k},i}(\widehat{\boldsymbol{\theta}}_{D})][\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\widehat{\boldsymbol{\theta}})+n^{-1}\sum_{k=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},k}(\widehat{\boldsymbol{\theta}})\psi_{u_{D,k},i}(\widehat{\boldsymbol{\theta}}_{D})]^{{\sf\scriptscriptstyle{T}}}$.

##  3 Numerical studies

###  3.1 Simulation studies

##### Data generation

We generated baseline covariates
$Z_{1,i}\sim\mathcal{N}(1,0.5)\times\mathbbm{1}[0,2]$, a truncated normal
distribution, and $Z_{2,i}\sim\mathcal{B}(0.8)$, a Bernoulli distribution.
Given ${\bf Z}_{i}=(Z_{1,i},Z_{2,i})^{{\sf\scriptscriptstyle{T}}}$, the
outcomes $(T,D)$ were generated from
$\log(T_{i}/3)=-\beta_{T,1}Z_{1,i}-\beta_{T,2}Z_{2,i}+e_{T,i}$ and
$\log(D_{i}/3)=-\beta_{D,1}Z_{1,i}-\beta_{D,2}Z_{2,i}+e_{D,i}$, where
$(e_{T,i},e_{D,i})$ followed a bivariate distribution that corresponds to a
certain parametric copula function. Here, we considered two commonly used
copula families: Gumbel and Clayton, with different strengths of association
specified by Kendall’s $\tau$, defined as
$\tau=4\int_{[0,1]^{2}}\mathbbm{C}(u_{1},u_{2})\mathrm{d}\mathbbm{C}(u_{1},u_{2})-1$.
Kendall’s $\tau$ has a one-to-one correspondence with $\alpha$: for example,
$\tau\equiv\tau(\alpha)=(\alpha-1)/\alpha$ for Gumbel and
$\tau(\alpha)=\alpha/(\alpha+2)$ for Clayton. We considered three dependence
levels $\tau=0.4,0.6,0.8$. As an example, to generate $(e_{T,i},e_{D,i})$ from
the Gumbel copula with $\tau=0.6$, we first generated $(u_{T,i},u_{D,i})$ with
parameter $\alpha=2.5$ corresponding to $\tau=0.6$, and then obtained
$e_{j,i}=\log(\log(u_{j,i}))$ for $j=T,D$. We set
$\boldsymbol{\beta}_{D}=(0.2,0)^{{\sf\scriptscriptstyle{T}}}$, which leads to
about 18% of death times censored, and considered two sets of
$\boldsymbol{\beta}_{T}=(1,1)^{{\sf\scriptscriptstyle{T}}},(1,0.5)^{{\sf\scriptscriptstyle{T}}}$,
which leads to different censoring rates, 3%-12%, and 5%-18%, receptively, for
the non-terminal event.

##### Simulation study I

The first simulation study focused on evaluating the estimation accuracy and
efficiency of the proposed PMLE approach and comparing the finite-sample
performances of the PMLE with the one-stage MLE method from Chen, (2012) under
each combination of $(\tau,\boldsymbol{\beta}_{T})$. The sample size $n$ was
set to 200 and 400, with 1,000 replications of simulated datasets.

We specified $G_{T}(t)=G_{D}(t)=t$ (i.e., PH model) for the marginal
distributions and the copula function as the true copula family. Both the PMLE
and one-stage MLE were obtained from maximizing the same log-likelihood
function (8) using the R package trust, which implements the trust region
optimization algorithm (Fletcher, 2000). For most combinations of the copula
family and $(n,\tau,\boldsymbol{\beta}_{T})$, the optimization converged for
all 1,000 replications. The exceptions are under Clayton with Kendall’s
$\tau=0.8$, there were seven and four non-converged replications for
$\boldsymbol{\beta}_{T}=(1,1)^{{\sf\scriptscriptstyle{T}}},(1,0.5)^{{\sf\scriptscriptstyle{T}}}$
respectively with $n=200$, and three and one non-converged replications for
$\boldsymbol{\beta}_{T}=(1,1)^{{\sf\scriptscriptstyle{T}}},(1,0.5)^{{\sf\scriptscriptstyle{T}}}$
respectively with $n=400$.

Tables 1 and 2 present the bias (BIAS), empirical standard deviation (ESD),
and the root mean square error (rMSE), relative to the true parameter value,
for the MLE and PMLE of $\alpha$ and $\boldsymbol{\beta}_{T}$ under Gumbel and
Clayton, respectively. For both copula families, the PMLE performs well with a
small bias, and the bias decreases as $n$ increases, demonstrating its
consistency. In addition, since a smaller sample size or a higher censoring
rate indicates less available information, the PMLE for the sample size
$n=200$ or $\boldsymbol{\beta}_{T}=(1,0.5)^{{\sf\scriptscriptstyle{T}}}$
(resulting in a higher censoring rate for $T$) produces larger BIAS, ESD, and
rMSE values for all the parameters. The same pattern is also observed for the
MLE.

Comparing the PMLE with MLE for the regression coefficients
$\boldsymbol{\beta}_{T}$, their biases are similar for $\beta_{T,1}$, but the
PMLE yields a smaller BIAS in estimating $\beta_{T,2}$. On the other hand, the
MLE is more efficient than the PMLE, as expected. For the copula parameter
$\alpha$, the PMLE has considerably smaller BIAS and ESD than the MLE.
Regarding the computational time, on average, for sample size $n=200$ and
$400$, respectively, the computational time of the MLE was 2.8 and 4.2 times
the time for the PMLE under Clayton, and 3.3 and 4.9 times under Gumbel. It
shows that the PMLE is computationally more efficient than the MLE while
maintaining comparable performance.

Table 3 presents the ratio of the ESD of the PMLE to its average standard
error (ASE) from the proposed analytic variance estimator. These ratios are
reasonably close to one, indicating the variance estimator is valid.

##### Simulation study II

We conducted a second simulation study to assess the robustness of the PMLE
and MLE under the misspecification of the copula function. Specifically, we
considered two settings: (i) we generated $(T,D)$ from the Gumbel copula but
fit the Clayton copula, and (ii) vice versa. The sample size was $n=400$,
Kendall’s $\tau=0.6$, and
$\boldsymbol{\beta}_{T}=(1,0.5)^{{\sf\scriptscriptstyle{T}}}$. When the copula
function is misspecified, the MLE or PMLE estimator $\widehat{\alpha}$ will
not converge to the true value of the copula parameter under the true copula
family, and the plug-in estimate of Kendall’s $\tau$,
$\widehat{\tau}=\tau(\widehat{\alpha})$, will not converge to the true
Kendall’s $\tau$ either. Thus, it is more meaningful to focus on the
estimation of the marginals in this setting. Table 4 presents the BIAS, ESD,
and rMSE for $\boldsymbol{\beta}_{D}$ and $\boldsymbol{\beta}_{T}$ with or
without the copula misspecification. When the copula is misspecified, the bias
for the PMLE of $\boldsymbol{\beta}_{D}$ is almost the same as the bias when
the copula is correctly specified. It is because our PMLE of
$\boldsymbol{\beta}_{D}$ is obtained without the information on the copula
function and, consequently, is still a consistent estimator under copula
misspecification. However, the bias for the MLE of $\boldsymbol{\beta}_{D}$
increases substantially under the copula misspecification. For the estimation
of $\boldsymbol{\beta}_{T}$, both the MLE and PMLE perform less ideally due to
the misspecification, but when the true copula is Gumbel but the Clayton
copula is fitted, the PMLE gives a much smaller bias, demonstrating higher
finite-sample robustness. When the true copula is Clayton but the Gumbel
copula is fitted, the biases for the PMLE and MLE are similar.

##### Simulation study III

In addition, we compared our PMLE with the two-stage estimator developed in
Zhu et al., (2021) under the Gumbel copula. Specifically, Zhu et al., (2021),
in the first stage, estimated the regression coefficient
$\boldsymbol{\beta}_{D}$ first, which produced
$\widehat{U}_{D,i}=S_{D}(X_{D,i}|{\bf
Z}_{i};\widehat{\boldsymbol{\beta}}_{D})$, and with the estimated marginal of
$D$ and a non-parametric estimator of the marginal of $T^{*}$ (time to the
first event, which is only subject to independent censoring), a PMLE of the
copula parameter was obtained. In the second stage, the time-varying
regression coefficient, $\boldsymbol{\beta}_{T}(t)$, was estimated via a
separate set of non-linear estimating equations. For comparison purposes, we
followed the same simulation setting of Zhu et al., (2021). Since the marginal
specifications are different, we focused on the comparison of the copula
parameter estimation in terms of the BIAS, ESD, and rMSE, presented in Table
5, where the results of Zhu et al., (2021)’s method were extracted from Table
1 of their paper directly. It needs to be pointed out that the definition of
the Gumbel copula parameter in Zhu et al., (2021) is different from ours, and
their Gumbel copula function is defined as

$\mathbbm{C}(u_{1},u_{2};\alpha^{*})=\exp\left(-\left[\\{-\log(u_{1})\\}^{1/\alpha^{*}}+\\{-\log(u_{2})\\}^{1/\alpha^{*}}\right]^{\alpha^{*}}\right),$

which indicates that our Gumbel copula parameter in (4) $\alpha=1/\alpha^{*}$.
Based on this relationship, we obtained the PMLE of $\alpha^{*}$ by
$\widehat{\alpha}^{*}=1/\widehat{\alpha}$, and with this estimator, we
calculated the summary statistics presented in Table 5. Across all the
combinations $(n,\tau,\boldsymbol{\beta}_{T})$, our PMLE has given more
desirable performances in all three metrics.

###  3.2 Analysis of the BMT data

We applied the proposed method to the BMT data, available in the R package
SemiCompRisks, which consists of 137 patients with acute myelocytic leukemia
(AML) or acute lymphoblastic leukemia (ALL) aged 7 to 52 (Klein and
Moeschberger, 2003). This data set contains time to death (terminal event,
$D$), which is subject to the administrative censoring of study time, and time
to relapse (non-terminal event, $T$), which is subject to censoring by both
death and study time. We were interested in investigating the dependence
between times to relapse and death, and the effects of disease groups, AML-
low, AML-high, and ALL, on the time to relapse and death. Here, AML-low was
set as the baseline group.

The PH models were assumed for both marginal distributions, and three copula
families, Clayton, Gumbel, and Frank, were considered. In addition, we assumed
that under each copula family, the copula parameter may vary with the disease
group. Both the MLE and PMLE approaches were implemented for estimating the
parameters. For inference, the standard errors of PMLE were calculated using
the proposed analytic variance estimators, and the standard errors of MLE were
calculated by a robust variance-covariance matrix estimator:
$\mathcal{I}_{n}^{-1}(\widehat{\boldsymbol{\theta}}_{\mbox{\tiny
MLE}})\mathcal{V}_{n}(\widehat{\boldsymbol{\theta}}_{\mbox{\tiny
MLE}})\mathcal{I}_{n}^{-1}(\widehat{\boldsymbol{\theta}}_{\mbox{\tiny MLE}})$
with
$\mathcal{I}_{n}(\boldsymbol{\theta})=-n^{-1}\sum_{i=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}\boldsymbol{\theta},i}(\boldsymbol{\theta})$
and
$\mathcal{V}_{n}(\boldsymbol{\theta})=n^{-1}\sum_{i=1}^{n}\dot{\ell}_{\boldsymbol{\theta},i}(\boldsymbol{\theta})\allowbreak\dot{\ell}_{\boldsymbol{\theta},i}(\boldsymbol{\theta})^{{\sf\scriptscriptstyle{T}}}$,
where
$\dot{\ell}_{\boldsymbol{\theta},i}(\boldsymbol{\theta})=\partial\ell_{i}(\boldsymbol{\theta})/\partial\boldsymbol{\theta}$
and
$\ddot{\ell}_{\boldsymbol{\theta}\boldsymbol{\theta},i}(\boldsymbol{\theta})=\partial^{2}\ell_{i}(\boldsymbol{\theta})/\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{{\sf\scriptscriptstyle{T}}}$
are the first-order and second-order derivatives of the log-likelihood
function $\ell_{i}(\boldsymbol{\theta})$ in (8) with respect to
$\boldsymbol{\theta}$.

The parameter estimates with their standard errors (in the parentheses) are
shown in Table 6. For the dependence between the bivariate event times, we
calculated the plug-in estimate of Kendall’s $\tau$:
$\widehat{\tau}=\tau(\widehat{\alpha})$ under each copula family, where
$\widehat{\alpha}$ is the MLE or PMLE of $\alpha$. As described earlier,
$\tau(\alpha)=(\alpha-1)/\alpha$ for Gumbel, and
$\tau(\alpha)=\alpha/(\alpha+2)$ for Clayton. In addition, for Frank,
$\tau(\alpha)=1-\frac{4}{\alpha}+\frac{4}{\alpha}\int_{0}^{\alpha}\frac{t}{e^{t}-1}\mathrm{d}t$.
The standard error of $\widehat{\tau}$ was obtained via the delta method:
$\widehat{\sigma}_{\tau}=|\dot{\tau}(\widehat{\alpha})|\widehat{\sigma}_{\alpha}$,
where $\widehat{\sigma}_{\tau}$ and $\widehat{\sigma}_{\alpha}$ denote the
standard error for $\widehat{\tau}$ and $\widehat{\alpha}$, respectively.

Comparing the MLE and PMLE, the two approaches produce similar point estimates
and standard errors. For marginal distributions, compared to the AML-low
group, both the AML-high and ALL groups have significantly higher risks of
death and relapse of leukemia under all the copula families. Figure 1 plots
the estimated survival function for the time to relapse for each disease group
with 95% confidence intervals. The AML-low group has the highest survival
rates, and the AML-high has the lowest rates. Regarding the dependence between
times to relapse and death, under each copula family, a strong and similar
association is observed for the three disease groups. A minor difference
across the three copula families is that under Clayton, the association is
strongest for the AML-low group, but the AML-high group has the strongest
association under Frank and Gumbel.

##  4 Conclusions

This article proposed a two-stage likelihood-based estimation procedure for
analyzing semi-competing risks data. The dependence between the non-terminal
and terminal event times is modeled by a parametric copula function, and the
marginals are specified via semi-parametric transformation models. One
motivation for such a stage-wise estimation is to alleviate the computational
challenge arising from the one-stage MLE approach, and through simulation
studies, our PMLE for the copula parameter has given a more desirable finite-
sample performance than the MLE. On the other hand, despite our method being
less efficient than the MLE for the regression coefficients of the marginals,
it is robust against copula misspecification. This is because the estimation
of $S_{D}(\cdot)$ in the first stage is guaranteed to be consistent, and as a
result, the estimation of $S_{T}(\cdot)$ might lead to a smaller bias than the
one-stage MLE.

Although our log-likelihood function in (8) is equivalent to the counting-
process-based log-likelihood in Chen, (2012) (Remark 3), the respective one-
stage MLEs would be different for finite samples due to the discretization of
$R_{j}(\cdot)$ for $j=T,D$. Thus, we observe close but not exactly the same
results between Table 6 of this article and Table 5 of Chen, (2012). However,
we expect this difference will decrease with a larger sample size.

The popular strategy for addressing the dependent censoring in estimating the
survival function of $T$ was using the first-event time $T^{*}$. However, its
implementation is limited to the Archimedean copula family, and when $T$ is
heavily censored by $D$, the estimation for the tail of $S_{T}(\cdot)$ is
challenging (Fine et al., 2001, Lakhal et al., 2008). In addition,
$S_{T^{*}}(\cdot)$ is often estimated non-parametrically even with the
covariates, which contributes to the less ideal estimation efficiency. Though
this could be improved by using a semi-parametric model, it is likely to
violate the model assumptions imposed on both $T$ and $D$. For example, if one
assumes that $T$ and $D$ both follow PH models, then the survival function of
$T^{*}$ does not follow a PH model and might not even have an explicit
expression.

Our Simulation Study II demonstrated that the estimation for $S_{T}(\cdot)$
and the dependence between $T$ and $D$ requires a correctly specified copula
function. The biggest challenge in designing a goodness-of-fit test is the
difficulty of distinguishing the misspecification of the copula function from
the misspecification of $S_{T}(\cdot)$, which is a future research problem.

## Acknowledgements

The work is supported by the grant DMS-2210481 from the National Science
Foundation.

## References

  * Chen, (2012) Chen, Y.-H. (2012).  Maximum likelihood analysis of semicompeting risks data with semiparametric regression models.  Lifetime Data Analysis, 18:36–57.
  * Fine et al., (2001) Fine, J. P., Jiang, H., and Chappell, R. (2001).  On semi-competing risks data.  Biometrika, 88(4):907–919.
  * Fletcher, (2000) Fletcher, R. (2000).  Practical Methods of Optimization.  John Wiley & Sons, Chichester, UK.
  * Ghosh, (2006) Ghosh, D. (2006).  Semiparametric inferences for association with semi-competing risks data.  Statistics in Medicine, 25(12):2059–2070.
  * Hsieh et al., (2008) Hsieh, J.-J., Wang, W., and Ding, A. A. (2008).  Regression analysis based on semicompeting risks data.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 70(1):3–20.
  * Klein and Moeschberger, (2003) Klein, J. P. and Moeschberger, M. L. (2003).  Survival Analysis: Techniques for Censored and Truncated Data.  Springer, New York, NY.
  * Lakhal et al., (2008) Lakhal, L., Rivest, L.-P., and Abdous, B. (2008).  Estimating survival and association in a semicompeting risks model.  Biometrics, 64(1):180–188.
  * Nelsen, (2006) Nelsen, R. B. (2006).  An Introduction to Copulas.  Springer, New York, NY, 2 edition.
  * Peng and Fine, (2007) Peng, L. and Fine, J. P. (2007).  Regression modeling of semicompeting risks data.  Biometrics, 63(1):96–108.
  * Shih and Louis, (1995) Shih, J. H. and Louis, T. A. (1995).  Inferences on the association parameter in copula models for bivariate survival data.  Biometrics, 51(4):1384–1399.
  * Sklar, (1959) Sklar, M. (1959).  Fonctions de répartition à n dimensions et leurs marges.  In Annales de l’ISUP, volume 8, pages 229–231.
  * van der Vaart and Wellner, (1996) van der Vaart, A. W. and Wellner, J. A. (1996).  Weak Convergence and Empirical Processes: With Applications to Statistics.  Springer, New York, NY.
  * Wang, (2003) Wang, W. (2003).  Estimating the association parameter for copula models under dependent censoring.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 65(1):257–273.
  * Zeng and Lin, (2006) Zeng, D. and Lin, D. (2006).  Efficient estimation of semiparametric transformation models for counting processes.  Biometrika, 93(3):627–640.
  * Zhou et al., (2016) Zhou, R., Zhu, H., Bondy, M., and Ning, J. (2016).  Semiparametric model for semi-competing risks data with application to breast cancer study.  Lifetime Data Analysis, 22:456–471.
  * Zhu et al., (2021) Zhu, H., Lan, Y., Ning, J., and Shen, Y. (2021).  Semiparametric copula-based regression modeling of semi-competing risks data.  Communications in Statistics-Theory and Methods, 51(22):7830–7845.

## Appendix

##  Appendix A Expressions of score functions

The score function $\Psi(\boldsymbol{\theta})$ is a
$(2p+\kappa_{T}+\kappa_{D}+1)$-dimensional vector, given by
$\Psi=\sum_{i=1}^{n}\Psi_{i}$, where

$\displaystyle\Psi_{i}=\frac{\partial}{\partial\boldsymbol{\theta}}\ell_{i}=\left(\begin{array}[]{ccccc}\Psi_{\alpha,i}&\Psi_{\boldsymbol{\beta}_{T},i}^{{\sf\scriptscriptstyle{T}}}&\Psi_{\mathrm{d}R_{T},i}^{{\sf\scriptscriptstyle{T}}}&\Psi_{\boldsymbol{\beta}_{D},i}^{{\sf\scriptscriptstyle{T}}}&\Psi_{\mathrm{d}R_{D},i}^{{\sf\scriptscriptstyle{T}}}\end{array}\right)^{{\sf\scriptscriptstyle{T}}}.$

The elements $\Psi_{\alpha,i}$ is scalar, $\Psi_{\boldsymbol{\beta}_{T},i}$
and $\Psi_{\boldsymbol{\beta}_{D},i}$ are $p$-dimensional vectors,
$\Psi_{\mathrm{d}R_{T},i}=(\Psi_{\mathrm{d}R_{T1},i},\ldots,\allowbreak\Psi_{\mathrm{d}R_{T\kappa_{T}},i})^{{\sf\scriptscriptstyle{T}}}$
is an $\kappa_{T}$-dimensional vector, and
$\Psi_{\mathrm{d}R_{D},i}=(\Psi_{\mathrm{d}R_{D1},i},\ldots,\Psi_{\mathrm{d}R_{D\kappa_{D}},i})^{{\sf\scriptscriptstyle{T}}}$
is a $\kappa_{D}$-dimensional vector.

To derive their expressions, we define the following notation: let
$\dot{f}(x)\equiv\mathrm{d}f(x)/\mathrm{d}x$ and
$\ddot{f}(x)\equiv\mathrm{d}^{2}f(x)/\mathrm{d}x^{2}$ denote the first-order
and second-order derivatives of a univariate function $f$. In addition,
$\dot{\mathbbm{C}}_{\alpha}=\partial\mathbbm{C}/\partial\alpha$,
$\ddot{\mathbbm{C}}_{u_{j}\alpha}=\partial^{2}\mathbbm{C}/\partial
u_{j}\partial\alpha$, and
$\ddot{\mathbbm{C}}_{u_{j}^{2}}=\partial^{2}\mathbbm{C}\partial u_{j}^{2}$ for
$j=T,D$, and
$\dddot{\mathbbm{C}}_{u_{T}^{\mbox{\tiny$d_{1}$}}u_{D}^{\mbox{\tiny$d_{2}$}}\alpha^{\mbox{\tiny$d_{3}$}}}=\partial^{3}\mathbbm{C}/\partial
u_{T}^{\mbox{\tiny$d_{1}$}}\partial
u_{D}^{\mbox{\tiny$d_{2}$}}\partial\alpha^{\mbox{\tiny$d_{3}$}}$ for integers
$d_{1},d_{2},d_{3}\geq 0$ with $d_{1}+d_{2}+d_{3}=3$. In the remainder of the
Appendix, for simplicity, we omit $(U_{T,i},U_{D,i};\alpha)$ from the copula
function $\mathbbm{C}$ as well as its partial derivatives.

The expression of $\Psi_{\alpha,i}$ is given as

$\Psi_{\alpha,i}=\delta_{T,i}\delta_{D,i}\frac{\dddot{\mathbbm{C}}_{u_{T}u_{D}\alpha}}{\ddot{\mathbbm{C}}_{u_{T}u_{D}}}+\delta_{T,i}(1-\delta_{D,i})\frac{\ddot{\mathbbm{C}}_{u_{T}\alpha}}{\dot{\mathbbm{C}}_{u_{T}}}+(1-\delta_{T,i})\delta_{D,i}\frac{\ddot{\mathbbm{C}}_{u_{D}\alpha}}{\dot{\mathbbm{C}}_{u_{D}}}+(1-\delta_{T,i})(1-\delta_{D,i})\frac{\dot{\mathbbm{C}}_{\alpha}}{\mathbbm{C}}.$

Define

$\displaystyle\Psi_{u_{T},i}$
$\displaystyle=\delta_{T,i}\delta_{D,i}\frac{\dddot{\mathbbm{C}}_{u_{T}^{2}u_{D}}}{\ddot{\mathbbm{C}}_{u_{T}u_{D}}}+\delta_{T,i}(1-\delta_{D,i})\frac{\ddot{\mathbbm{C}}_{u_{T}^{2}}}{\dot{\mathbbm{C}}_{u_{T}}}+(1-\delta_{T,i})\delta_{D,i}\frac{\ddot{\mathbbm{C}}_{u_{T}u_{D}}}{\dot{\mathbbm{C}}_{u_{D}}}+(1-\delta_{T,i})(1-\delta_{D,i})\frac{\dot{\mathbbm{C}}_{u_{T}}}{\mathbbm{C}},$

and

$\displaystyle\Psi_{u_{D},i}$
$\displaystyle=\delta_{T,i}\delta_{D,i}\frac{\dddot{\mathbbm{C}}_{u_{T}u_{D}^{2}}}{\ddot{\mathbbm{C}}_{u_{T}u_{D}}}+\delta_{T,i}(1-\delta_{D,i})\frac{\ddot{\mathbbm{C}}_{u_{T}u_{D}}}{\dot{\mathbbm{C}}_{u_{T}}}+(1-\delta_{T,i})\delta_{D,i}\frac{\ddot{\mathbbm{C}}_{u_{D}^{2}}}{\dot{\mathbbm{C}}_{u_{D}}}+(1-\delta_{T,i})(1-\delta_{D,i})\frac{\dot{\mathbbm{C}}_{u_{D}}}{\mathbbm{C}}.$

Let $R_{T,i}=R_{T}(X_{i})$ and
$\Lambda_{T,i}=R_{T,i}e^{\boldsymbol{\beta}_{T}^{\prime}{\bf Z}_{i}}$. In
addition, let $R_{D,i}=R_{D}(C_{i})$ and
$\Lambda_{iD}=R_{D,i}e^{\boldsymbol{\beta}_{D}^{\prime}{\bf Z}_{i}}$. The
expressions of $\Psi_{\boldsymbol{\beta}_{T},i}$ and
$\Psi_{\boldsymbol{\beta}_{D},i}$ are, for $j=T,D$,

$\Psi_{\boldsymbol{\beta}_{j},i}={\bf
Z}_{i}\left\\{-\Psi_{u_{j},i}U_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)R_{j,i}e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}-\delta_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)R_{j,i}e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}+\delta_{j,i}\frac{\ddot{G}_{j}\left(\Lambda_{j,i}\right)}{\dot{G}_{j}\left(\Lambda_{j,i}\right)}R_{j,i}e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}+\delta_{j,i}\right\\}.$

Let $Y_{T,l,i}=\mathbbm{1}(X_{i}\geq t_{l})$ and
$I_{T,l,i}=\mathbbm{1}(X_{i}=t_{l})$. In addition, let
$Y_{D,l,i}=\mathbbm{1}(C_{i}\geq d_{l})$ and
$I_{D,l,i}=\mathbbm{1}(C_{i}=d_{l})$. The expressions of
$\Psi_{\mathrm{d}R_{T,l},i}$ and $\Psi_{\mathrm{d}R_{D,l},i}$ follow a similar
form. For $j=T,D$,

$\Psi_{\mathrm{d}R_{j,l},i}=Y_{j,l,i}\left\\{-\Psi_{u_{j},i}U_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}-\delta_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}+\delta_{j,i}\frac{\ddot{G}_{j}\left(\Lambda_{j,i}\right)}{\dot{G}_{j}\left(\Lambda_{j,i}\right)}e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\right\\}+\delta_{j,i}I_{j,l,i}\mathrm{d}R_{j,l}^{-1}.$

##  Appendix B Expressions of information matrices

The full information matrix
$\mathcal{I}_{n}(\boldsymbol{\theta})=n^{-1}\sum_{i=1}^{n}\mathcal{I}_{i}(\boldsymbol{\theta})$
is a $(2p+\kappa_{D}+\kappa_{T}+1)$-dimensional square matrix, where

$\displaystyle\mathcal{I}_{i}(\boldsymbol{\theta})=-\frac{\partial^{2}}{\partial\boldsymbol{\theta}\partial\boldsymbol{\theta}^{{\sf\scriptscriptstyle{T}}}}\ell_{i}(\boldsymbol{\theta})=\left[\begin{array}[]{ccccc}\mathcal{I}_{\alpha\alpha,i}&\mathcal{I}_{\alpha\boldsymbol{\beta}_{T},i}^{{\sf\scriptscriptstyle{T}}}&\mathcal{I}_{\alpha\mathrm{d}R_{T},i}^{{\sf\scriptscriptstyle{T}}}&\mathcal{I}_{\alpha\boldsymbol{\beta}_{D},i}^{{\sf\scriptscriptstyle{T}}}&\mathcal{I}_{\alpha\mathrm{d}R_{D},i}^{{\sf\scriptscriptstyle{T}}}\\\
\mathcal{I}_{\alpha\boldsymbol{\beta}_{T},i}&\mathcal{I}_{\boldsymbol{\beta}_{T}\boldsymbol{\beta}_{T},i}&\mathcal{I}_{\boldsymbol{\beta}_{T}\mathrm{d}R_{T},i}^{{\sf\scriptscriptstyle{T}}}&\mathcal{I}_{\boldsymbol{\beta}_{T}\boldsymbol{\beta}_{D},i}^{{\sf\scriptscriptstyle{T}}}&\mathcal{I}_{\boldsymbol{\beta}_{T}\mathrm{d}R_{D},i}^{{\sf\scriptscriptstyle{T}}}\\\
\mathcal{I}_{\alpha\mathrm{d}R_{T},i}&\mathcal{I}_{\mathrm{d}R_{T}\boldsymbol{\beta}_{T},i}&\mathcal{I}_{\mathrm{d}R_{T}\mathrm{d}R_{T},i}&\mathcal{I}_{\mathrm{d}R_{T}\boldsymbol{\beta}_{D},i}^{{\sf\scriptscriptstyle{T}}}&\mathcal{I}_{\mathrm{d}R_{T}\mathrm{d}R_{D},i}^{{\sf\scriptscriptstyle{T}}}\\\
\mathcal{I}_{\alpha\boldsymbol{\beta}_{D}}&\mathcal{I}_{\boldsymbol{\beta}_{D}\boldsymbol{\beta}_{T},i}&\mathcal{I}_{\boldsymbol{\beta}_{D}\mathrm{d}R_{T},i}&\mathcal{I}_{\boldsymbol{\beta}_{D}\boldsymbol{\beta}_{D},i}&\mathcal{I}_{\boldsymbol{\beta}_{D}\mathrm{d}R_{D},i}^{{\sf\scriptscriptstyle{T}}}\\\
\mathcal{I}_{\alpha\mathrm{d}R_{D}}&\mathcal{I}_{\mathrm{d}R_{D}\boldsymbol{\beta}_{T},i}&\mathcal{I}_{\mathrm{d}R_{D}\mathrm{d}R_{T},i}&\mathcal{I}_{\mathrm{d}R_{D}\boldsymbol{\beta}_{D},i}&\mathcal{I}_{\mathrm{d}R_{D}\mathrm{d}R_{D},i}\end{array}\right].$

To derive the expression of the elements in
$\mathcal{I}_{i}(\boldsymbol{\theta})$, we need to define let
$\ddddot{\mathbbm{C}}_{u_{T}^{\mbox{\tiny$d_{1}$}}u_{D}^{\mbox{\tiny$d_{2}$}}\alpha^{\mbox{\tiny$d_{3}$}}}=\partial^{4}\mathbbm{C}/\partial
u_{T}^{\mbox{\tiny$d_{1}$}}\partial
u_{D}^{\mbox{\tiny$d_{2}$}}\partial\alpha^{\mbox{\tiny$d_{3}$}}$ for integers
$d_{1},d_{2},d_{3}\geq 0$ with $d_{1}+d_{2}+d_{3}=4$. We have
$\mathcal{I}_{\alpha\alpha,i}$, a scalar, given by

$\displaystyle\mathcal{I}_{\alpha\alpha,i}$
$\displaystyle=\delta_{T,i}\delta_{D,i}\frac{\ddddot{\mathbbm{C}}_{u_{T}u_{D}\alpha^{2}}\ddot{\mathbbm{C}}_{u_{T}u_{D}}-\dddot{\mathbbm{C}}_{u_{T}u_{D}\alpha}^{2}}{\ddot{\mathbbm{C}}_{u_{T}u_{D}}^{2}}+\delta_{T,i}(1-\delta_{D,i})\frac{\dddot{\mathbbm{C}}_{u_{T}\alpha^{2}}\dot{\mathbbm{C}}_{u_{T}}-\ddot{\mathbbm{C}}_{u_{T}\alpha}^{2}}{\dot{\mathbbm{C}}_{u_{T}}^{2}}+$
$\displaystyle\quad(1-\delta_{T,i})\delta_{D,i}\frac{\dddot{\mathbbm{C}}_{u_{D}\alpha^{2}}\dot{\mathbbm{C}}_{u_{D}}-\ddot{\mathbbm{C}}_{u_{D}\alpha}^{2}}{\dot{\mathbbm{C}}_{u_{D}}^{2}}+(1-\delta_{T,i})(1-\delta_{D,i})\frac{\ddot{\mathbbm{C}}_{\alpha^{2}}\mathbbm{C}-\dot{\mathbbm{C}}_{\alpha}^{2}}{\mathbbm{C}^{2}}.$

Define the following

$\displaystyle\mathcal{I}_{u_{T}\alpha,i}$
$\displaystyle=\delta_{T,i}\delta_{D,i}\frac{\ddddot{\mathbbm{C}}_{u_{T}^{2}u_{D}\alpha}\ddot{\mathbbm{C}}_{u_{T}u_{D}}-\dddot{\mathbbm{C}}_{u_{T}^{2}u_{D}}\dddot{\mathbbm{C}}_{u_{T}u_{D}\alpha}}{\ddot{\mathbbm{C}}_{u_{T}u_{D}}^{2}}+\delta_{T,i}(1-\delta_{D,i})\frac{\dddot{\mathbbm{C}}_{u_{T}^{2}\alpha}\dot{\mathbbm{C}}_{u_{T}}-\ddot{\mathbbm{C}}_{u_{T}^{2}}\ddot{\mathbbm{C}}_{u_{T}\alpha}}{\dot{\mathbbm{C}}_{u_{T}}^{2}}+$
$\displaystyle\quad(1-\delta_{T,i})\delta_{D,i}\frac{\dddot{\mathbbm{C}}_{u_{T}u_{D}\alpha}\dot{\mathbbm{C}}_{u_{D}}-\ddot{\mathbbm{C}}_{u_{T}u_{D}}\ddot{\mathbbm{C}}_{u_{D}\alpha}}{\dot{\mathbbm{C}}_{u_{D}}^{2}}+(1-\delta_{T,i})(1-\delta_{D,i})\frac{\ddot{\mathbbm{C}}_{u_{T}\alpha}\mathbbm{C}-\dot{\mathbbm{C}}_{u_{T}}\dot{\mathbbm{C}}_{\alpha}}{\mathbbm{C}^{2}},$
$\displaystyle\mathcal{I}_{u_{D}\alpha,i}$
$\displaystyle=\delta_{T,i}\delta_{D,i}\frac{\ddddot{\mathbbm{C}}_{u_{T}u_{D}^{2}\alpha}\ddot{\mathbbm{C}}_{u_{T}u_{D}}-\dddot{\mathbbm{C}}_{u_{T}u_{D}^{2}}\dddot{\mathbbm{C}}_{u_{T}u_{D}\alpha}}{\ddot{\mathbbm{C}}_{u_{T}u_{D}}^{2}}+\delta_{T,i}(1-\delta_{D,i})\frac{\dddot{\mathbbm{C}}_{u_{T}u_{D}\alpha}\dot{\mathbbm{C}}_{u_{T}}-\ddot{\mathbbm{C}}_{u_{T}u_{D}}\ddot{\mathbbm{C}}_{u_{T}\alpha}}{\dot{\mathbbm{C}}_{u_{T}}^{2}}+$
$\displaystyle\quad(1-\delta_{T,i})\delta_{D,i}\frac{\dddot{\mathbbm{C}}_{u_{D}^{2}\alpha}\dot{\mathbbm{C}}_{u_{D}}-\ddot{\mathbbm{C}}_{u_{D}^{2}}\ddot{\mathbbm{C}}_{u_{D}\alpha}}{\dot{\mathbbm{C}}_{u_{D}}^{2}}+(1-\delta_{T,i})(1-\delta_{D,i})\frac{\ddot{\mathbbm{C}}_{u_{D}\alpha}\mathbbm{C}-\dot{\mathbbm{C}}_{u_{D}}\dot{\mathbbm{C}}_{\alpha}}{\mathbbm{C}^{2}},$
$\displaystyle\mathcal{I}_{u_{T}u_{T},i}$
$\displaystyle=\delta_{T,i}\delta_{D,i}\frac{\ddddot{\mathbbm{C}}_{u_{T}^{3}u_{D}}\ddot{\mathbbm{C}}_{u_{T}u_{D}}-\dddot{\mathbbm{C}}_{u_{T}^{2}u_{D}}^{2}}{\ddot{\mathbbm{C}}_{u_{T}u_{D}}^{2}}+\delta_{T,i}(1-\delta_{D,i})\frac{\dddot{\mathbbm{C}}_{u_{T}^{3}}\dot{\mathbbm{C}}_{u_{T}}-\ddot{\mathbbm{C}}_{u_{T}^{2}}^{2}}{\dot{\mathbbm{C}}_{u_{T}}^{2}}+$
$\displaystyle\quad(1-\delta_{T,i})\delta_{D,i}\frac{\dddot{\mathbbm{C}}_{u_{T}^{2}u_{D}}\dot{\mathbbm{C}}_{u_{D}}-\ddot{\mathbbm{C}}_{u_{T}u_{D}}^{2}}{\dot{\mathbbm{C}}_{u_{D}}^{2}}+(1-\delta_{T,i})(1-\delta_{D,i})\frac{\ddot{\mathbbm{C}}_{u_{T}^{2}}\mathbbm{C}-\dot{\mathbbm{C}}_{u_{T}}^{2}}{\mathbbm{C}^{2}},$
$\displaystyle\mathcal{I}_{u_{D}u_{D},i}$
$\displaystyle=\delta_{T,i}\delta_{D,i}\frac{\ddddot{\mathbbm{C}}_{u_{T}u_{D}^{3}}\ddot{\mathbbm{C}}_{u_{T}u_{D}}-\dddot{\mathbbm{C}}_{u_{T}u_{D}^{2}}^{2}}{\ddot{\mathbbm{C}}_{u_{T}u_{D}}^{2}}+\delta_{T,i}(1-\delta_{D,i})\frac{\dddot{\mathbbm{C}}_{u_{T}u_{D}^{2}}\dot{\mathbbm{C}}_{u_{T}}-\ddot{\mathbbm{C}}_{u_{T}u_{D}}^{2}}{\dot{\mathbbm{C}}_{u_{T}}^{2}}+$
$\displaystyle\quad(1-\delta_{T,i})\delta_{D,i}\frac{\dddot{\mathbbm{C}}_{u_{D}^{3}}\dot{\mathbbm{C}}_{u_{D}}-\ddot{\mathbbm{C}}_{u_{D}^{2}}^{2}}{\dot{\mathbbm{C}}_{u_{D}}^{2}}+(1-\delta_{T,i})(1-\delta_{D,i})\frac{\ddot{\mathbbm{C}}_{u_{D}^{2}}\mathbbm{C}-\dot{\mathbbm{C}}_{u_{D}}^{2}}{\mathbbm{C}^{2}},$

and

$\displaystyle\mathcal{I}_{u_{T}u_{D},i}$
$\displaystyle=\delta_{T,i}\delta_{D,i}\frac{\ddddot{\mathbbm{C}}_{u_{T}^{2}u_{D}^{2}}\ddot{\mathbbm{C}}_{u_{T}u_{D}}-\dddot{\mathbbm{C}}_{u_{T}u_{D}^{2}}\dddot{\mathbbm{C}}_{u_{T}^{2}u_{D}}}{\ddot{\mathbbm{C}}_{u_{T}u_{D}}^{2}}+\delta_{T,i}(1-\delta_{D,i})\frac{\dddot{\mathbbm{C}}_{u_{T}^{2}u_{D}}\dot{\mathbbm{C}}_{u_{T}}-\ddot{\mathbbm{C}}_{u_{T}u_{D}}\ddot{\mathbbm{C}}_{u_{T}^{2}}}{\dot{\mathbbm{C}}_{u_{T}}^{2}}+$
$\displaystyle\quad(1-\delta_{T,i})\delta_{D,i}\frac{\dddot{\mathbbm{C}}_{u_{T}u_{D}^{2}}\dot{\mathbbm{C}}_{u_{D}}-\ddot{\mathbbm{C}}_{u_{D}^{2}}\ddot{\mathbbm{C}}_{u_{T}u_{D}}}{\dot{\mathbbm{C}}_{u_{D}}^{2}}+(1-\delta_{T,i})(1-\delta_{D,i})\frac{\ddot{\mathbbm{C}}_{u_{T}u_{D}}\mathbbm{C}-\dot{\mathbbm{C}}_{u_{D}}\dot{\mathbbm{C}}_{u_{T}}}{\mathbbm{C}^{2}}.$

Both $\mathcal{I}_{\boldsymbol{\beta}_{T}\boldsymbol{\beta}_{T},i}$ and
$\mathcal{I}_{\boldsymbol{\beta}_{D}\boldsymbol{\beta}_{D},i}$ are
$p$-dimensional square matrices, and their expression follow similar forms.
For $j=T,D$,

$\displaystyle\mathcal{I}_{\boldsymbol{\beta}_{j}\boldsymbol{\beta}_{j},i}=$
$\displaystyle{\bf Z}_{i}{\bf
Z}_{i}^{{\sf\scriptscriptstyle{T}}}\Bigg{\\{}\mathcal{I}_{u_{j}u_{j},i}U_{j,i}^{2}\dot{G}_{j}\left(\Lambda_{j,i}\right)^{2}R_{j,i}^{2}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}+\Psi_{u_{j},i}U_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)^{2}R_{j,i}^{2}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}$
$\displaystyle\qquad\quad-\Psi_{u_{j},i}U_{j,i}\ddot{G}_{j}\left(\Lambda_{j,i}\right)R_{j,i}^{2}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}-\Psi_{u_{j},i}U_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)R_{j,i}e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}$
$\displaystyle\qquad\quad-\delta_{j,i}\ddot{G}_{j}\left(\Lambda_{j,i}\right)R_{j,i}^{2}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}-\delta_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)R_{j,i}e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}$
$\displaystyle\qquad\quad+\delta_{j,i}\frac{\dddot{G}_{j}\left(\Lambda_{j,i}\right)\dot{G}_{j}\left(\Lambda_{j,i}\right)-\ddot{G}_{j}\left(\Lambda_{j,i}\right)^{2}}{\dot{G}_{j}\left(\Lambda_{j,i}\right)}R_{j,i}^{2}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}+\delta_{j,i}\frac{\ddot{G}_{j}\left(\Lambda_{j,i}\right)}{\dot{G}_{j}\left(\Lambda_{j,i}\right)}R_{j,i}e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\Bigg{\\}}.$

The expressions of $\mathcal{I}_{\mathrm{d}R_{T}\mathrm{d}R_{T},i}$, a
$\kappa_{T}\times\kappa_{T}$-dimensional matrix, and
$\mathcal{I}_{\mathrm{d}R_{D}\mathrm{d}R_{D},i}$, a
$\kappa_{D}\times\kappa_{D}$-dimensional matrix, follow similar forms. For
$j=T,D$ and $l,l^{\prime}=1,\cdots,\kappa_{j}$, the $(l,l^{\prime})$ element
of the matrix $\mathcal{I}_{\mathrm{d}R_{j}\mathrm{d}R_{j},i}$ is given by

$\displaystyle\mathcal{I}_{\mathrm{d}R_{j,l}\mathrm{d}R_{j,l^{\prime}},i}$
$\displaystyle=$ $\displaystyle
Y_{j,l,i}Y_{j,l^{\prime},i}\Bigg{\\{}\mathcal{I}_{u_{j}u_{j},i}U_{j,i}^{2}\dot{G}_{j}\left(\Lambda_{j,i}\right)^{2}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}+\Psi_{u_{j},i}U_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)^{2}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}-\Psi_{u_{j},i}U_{j,i}\ddot{G}_{j}\left(\Lambda_{j,i}\right)e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}$
$\displaystyle-\delta_{j,i}\ddot{G}_{j}\left(\Lambda_{j,i}\right)e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}+\delta_{j,i}\frac{\dddot{G}_{j}\left(\Lambda_{j,i}\right)\dot{G}_{j}\left(\Lambda_{j,i}\right)-\ddot{G}_{j}\left(\Lambda_{j,i}\right)^{2}}{\dot{G}_{j}\left(\Lambda_{j,i}\right)}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}-\delta_{j,i}I_{j,l,i}\mathbbm{1}(l=l^{\prime})\mathrm{d}R_{j,l}^{-2}\Bigg{\\}}.$

Both $\mathcal{I}_{\alpha\boldsymbol{\beta}_{T},i}$ and
$\mathcal{I}_{\alpha\boldsymbol{\beta}_{D},i}$ are $p$-dimensional vectors,
and their expressions follow similar forms. For $j=T,D$,

$\mathcal{I}_{\alpha\boldsymbol{\beta}_{j},i}=-\mathcal{I}_{u_{j}\alpha,i}U_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)R_{j,i}e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}{\bf Z}_{i}.$

The expressions of $\mathcal{I}_{\alpha\mathrm{d}R_{T},i}$, a
$\kappa_{T}$-dimensional vector, and $\mathcal{I}_{\alpha\mathrm{d}R_{D},i}$,
a $\kappa_{D}$-dimensional vector, follow similar forms. For $j=T,D$ and
$l=1,\ldots,\kappa_{j}$, the $l$-th element of the vector
$\mathcal{I}_{\alpha\mathrm{d}R_{j},i}$ is given as

$\mathcal{I}_{\alpha\mathrm{d}R_{j,l},i}=-Y_{j,l,i}\mathcal{I}_{u_{j}\alpha,i}U_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}.$

$\mathcal{I}_{\boldsymbol{\beta}_{T}\boldsymbol{\beta}_{D},i}$ is a $p\times
p$-dimensional matrix, given by

$\mathcal{I}_{\boldsymbol{\beta}_{T}\boldsymbol{\beta}_{D},i}={\bf Z}_{i}{\bf
Z}_{i}^{{\sf\scriptscriptstyle{T}}}\mathcal{I}_{u_{T}u_{D},i}U_{T,i}U_{D,i}e^{(\boldsymbol{\beta}_{T}+\boldsymbol{\beta}_{D})^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}R_{T,i}R_{D,i}\dot{G}_{T}\left(\Lambda_{T,i}\right)\dot{G}_{D}\left(\Lambda_{D,i}\right).$

$\mathcal{I}_{\mathrm{d}R_{T}\mathrm{d}R_{D},i}$ is a
$\kappa_{T}\times\kappa_{D}$-dimensional matrix with the $(l,l^{\prime})$
element, given by

$\mathcal{I}_{\mathrm{d}R_{T,l}\mathrm{d}R_{D,l^{\prime}},i}=Y_{T,l,i}Y_{D,l^{\prime},i}\mathcal{I}_{u_{T}u_{D},i}U_{T,i}U_{D,i}e^{(\boldsymbol{\beta}_{T}+\boldsymbol{\beta}_{D})^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\dot{G}_{T}\left(\Lambda_{T,i}\right)\dot{G}_{D}\left(\Lambda_{D,i}\right).$

The expressions of $\mathcal{I}_{\boldsymbol{\beta}_{T}\mathrm{d}R_{T},i}$, a
$p\times\kappa_{T}$-dimensional matrix, and
$\mathcal{I}_{\boldsymbol{\beta}_{D}\mathrm{d}R_{D},i}$, a
$p\times\kappa_{D}$-dimensional matrix, follow a similar form. For $j=T,D$ and
$l=1,\cdots\kappa_{j}$, the $l$-th column of the matrix
$\mathcal{I}_{\boldsymbol{\beta}_{j}\mathrm{d}R_{j},i}$ is given by

$\displaystyle\mathcal{I}_{\boldsymbol{\beta}_{j}\mathrm{d}R_{j,l},i}=$
$\displaystyle{\bf
Z}_{i}Y_{j,l,i}\Bigg{\\{}\mathcal{I}_{u_{j}u_{j},i}U_{j,i}^{2}\dot{G}_{j}\left(\Lambda_{j,i}\right)^{2}R_{j,i}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}+\Psi_{u_{j},i}U_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)^{2}R_{j,i}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}$
$\displaystyle\qquad\qquad-\Psi_{u_{j},i}U_{j,i}\ddot{G}_{j}\left(\Lambda_{j,i}\right)R_{j,i}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}-\Psi_{u_{j},i}U_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}$
$\displaystyle\qquad\qquad-\delta_{j,i}\ddot{G}_{j}\left(\Lambda_{j,i}\right)R_{j,i}^{2}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}-\delta_{j,i}\dot{G}_{j}\left(\Lambda_{j,i}\right)e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}$
$\displaystyle\qquad\qquad+\delta_{j,i}\frac{\dddot{G}_{j}\left(\Lambda_{j,i}\right)\dot{G}_{j}\left(\Lambda_{j,i}\right)-\ddot{G}_{j}\left(\Lambda_{j,i}\right)^{2}}{\dot{G}_{j}\left(\Lambda_{j,i}\right)}R_{j,i}e^{2\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}+\delta_{j,i}\frac{\ddot{G}_{j}\left(\Lambda_{j,i}\right)}{\dot{G}_{j}\left(\Lambda_{j,i}\right)}e^{\boldsymbol{\beta}_{j}^{{\sf\scriptscriptstyle{T}}}{\bf
Z}_{i}}\Bigg{\\}}.$

The expressions of $\mathcal{I}_{\boldsymbol{\beta}_{T}\mathrm{d}R_{D},i}$, a
$p\times\kappa_{D}$-dimensional matrix, and
$\mathcal{I}_{\boldsymbol{\beta}_{D}\mathrm{d}R_{T},i}$, a
$p\times\kappa_{T}$-dimensional matrix, follow similar forms. For
$j,j^{\prime}=T,D$ with $j\neq j^{\prime}$ and
$l=1,\cdots\kappa_{j^{\prime}}$, the $l$-th column of the matrix
$\mathcal{I}_{\boldsymbol{\beta}_{j}\mathrm{d}R_{j^{\prime}},i}$ is given by

$\mathcal{I}_{\boldsymbol{\beta}_{j}\mathrm{d}R_{j^{\prime},k},i}={\bf
Z}_{i}Y_{j^{\prime},k,i}U_{T,i}U_{D,i}e^{(\boldsymbol{\beta}_{T}+\boldsymbol{\beta}_{D})^{\prime}{\bf
Z}_{i}}\dot{G}_{T}(\Lambda_{T,i})\dot{G}_{T}(\Lambda_{D,i})R_{j,i}.$

##  Appendix C Proof of Theorems 1 and 2

The proof of Theorems 1 and 2 requires the asymptotic properties of
$\widehat{\boldsymbol{\theta}}_{D}$ obtained in the first stage, which was
established in Zeng and Lin, (2006), where it was shown that
$\widehat{\boldsymbol{\beta}}_{D}$ converges to $\boldsymbol{\beta}_{D}^{0}$
and $\widehat{R}_{D}(\cdot)$ converges to $R_{D}^{0}(\cdot)$ uniformly in the
interval $[0,\xi]$. In addition, by Theorem 3.3.1 of van der Vaart and
Wellner, (1996),
$\sqrt{n}(\widehat{\boldsymbol{\theta}}_{D}-\boldsymbol{\theta}_{D}^{0})$ can
be expressed as
$\sqrt{n}(\widehat{\boldsymbol{\theta}}_{D}-\boldsymbol{\theta}_{D}^{0})=n^{-1/2}\sum_{i=1}^{n}\psi_{\boldsymbol{\theta}_{D},i}(\boldsymbol{\theta}_{D}^{0})+o_{p}(1)$
with
$\psi_{\boldsymbol{\theta}_{D},i}(\boldsymbol{\theta}_{D}^{0})=\left(\mathcal{I}_{D}^{0}\right)^{-1}\Psi_{D,i}(\boldsymbol{\theta}_{D}^{0})$,
where
$\mathcal{I}_{D}^{0}=\mathbb{E}\\{\mathcal{I}_{D,i}(\boldsymbol{\theta}_{D}^{0})\\}$.
Since the variance-covariance matrix of
$n^{-1/2}\sum_{i=1}^{n}\Psi_{D,i}(\boldsymbol{\theta}_{D}^{0})$ is equal to
$\mathcal{I}_{D}^{0}$, by the functional delta method,
$\sqrt{n}\\{F_{D}(\widehat{\boldsymbol{\theta}}_{D})-F_{D}(\boldsymbol{\theta}_{D}^{0})\\}$
converges to a zero-mean Gaussian process with the variance-covariance matrix
$\dot{F}_{D}(\boldsymbol{\theta}_{D}^{0})^{{\sf\scriptscriptstyle{T}}}(\mathcal{I}_{D}^{0})^{-1}\dot{F}_{D}(\boldsymbol{\theta}_{D}^{0})$.

Define ${\mathfrak{H}}(\boldsymbol{\theta}_{D};t,{\bf
z})=\exp[-G_{D}\\{R_{D}(t)e^{\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
z}}\\}]$ and let

${\mathfrak{h}}(\boldsymbol{\theta}_{D};t,{\bf
z})=\frac{\partial{\mathfrak{H}}(\boldsymbol{\theta}_{D};t,{\bf
z})}{\partial\boldsymbol{\theta}_{D}}=\left({\mathfrak{h}}_{\boldsymbol{\beta}_{D}}^{{\sf\scriptscriptstyle{T}}},{\mathfrak{h}}_{\mathrm{d}R_{D,1}},\ldots,{\mathfrak{h}}_{\mathrm{d}R_{D,\kappa_{D}}}\right)^{{\sf\scriptscriptstyle{T}}},$

where

$\displaystyle{\mathfrak{h}}_{\boldsymbol{\beta}_{D}}$
$\displaystyle=-{\mathfrak{h}}(\boldsymbol{\theta}_{D};t,{\bf
z})\dot{G}_{D}\left\\{R_{D}(t)e^{\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
z}}\right\\}R_{D}(t){\bf z},$
$\displaystyle{\mathfrak{h}}_{\mathrm{d}R_{D,l}}$
$\displaystyle=-{\mathfrak{h}}(\boldsymbol{\theta}_{D};t,{\bf
z})\dot{G}_{D}\left\\{R_{D}(t)e^{\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
z}}\right\\}e^{\boldsymbol{\beta}_{D}^{{\sf\scriptscriptstyle{T}}}{\bf
z}}\mathbbm{1}(t\geq d_{l}),\,l=1,\ldots,\kappa_{D}.$

Let
$\widehat{U}_{D,i}={\mathfrak{H}}(\widehat{\boldsymbol{\theta}}_{D};C_{i},{\bf
Z}_{i})$ and
$U_{D,i}^{0}={\mathfrak{H}}(\boldsymbol{\theta}_{D}^{0};C_{i},{\bf Z}_{i})$.
By the functional delta method, we have
$\sqrt{n}(\widehat{U}_{D,i}-U_{D,i}^{0})={\mathfrak{h}}(\boldsymbol{\theta}_{D}^{0};C_{i},{\bf
Z}_{i})^{{\sf\scriptscriptstyle{T}}}\\{\sqrt{n}(\widehat{\boldsymbol{\theta}}_{D}-\boldsymbol{\theta}_{D}^{0})\\}+o_{p}(1)$.
And, by the expansion of
$\sqrt{n}(\widehat{\boldsymbol{\theta}}_{D}-\boldsymbol{\theta}_{D}^{0})$
above, we have

$\sqrt{n}(\widehat{U}_{D,i}-U_{D,i}^{0})=n^{-1}\sum_{k=1}^{n}\psi_{u_{D,i},k}(\boldsymbol{\theta}_{D}^{0})+o_{p}(1),$
(9)

where
$\psi_{u_{D,i},k}(\boldsymbol{\theta}_{D}^{0})={\mathfrak{h}}(\boldsymbol{\theta}_{D}^{0};C_{i},{\bf
Z}_{i})\psi_{\boldsymbol{\theta}_{D},k}(\boldsymbol{\theta}_{D}^{0})$.

###  C.1 Proof of Theorem 1

Let
$\breve{\boldsymbol{\theta}}_{1}=\arg\max_{\boldsymbol{\theta}_{1}}\ell(\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{D}^{0})$.
By similar arguments in Chen, (2012), we have $\breve{\alpha}$ converges
$\alpha^{0}$ and $\breve{\boldsymbol{\beta}}_{T}$ converges
$\boldsymbol{\beta}_{T}^{0}$, and
$\breve{R}_{T}(\cdot)\equiv\\{\breve{\mathrm{d}R}_{T,l},l=1,\ldots,\kappa_{T}\\}$
converges to $R_{T}^{0}(\cdot)$ uniformly in the interval $[0,\xi]$. By C4 and
the consistency of $\widehat{\boldsymbol{\theta}}_{D}$, we have, almost
surely,
$|\ell(\boldsymbol{\theta}_{1},\widehat{\boldsymbol{\theta}}_{D})-\ell(\boldsymbol{\theta}_{1},\boldsymbol{\theta}_{D}^{0})|_{\infty}=o_{p}(1)$.
Thus, $|\widehat{\alpha}-\breve{\alpha}|=o_{p}(1)$,
$\|\widehat{\boldsymbol{\beta}}_{T}-\breve{\boldsymbol{\beta}}_{T}\|=o_{p}(1)$,
and $|\widehat{R}_{T}-\breve{R}_{T}|_{\infty}=o_{p}(1)$. With probability 1,
$\widehat{\alpha}$ converges $\alpha^{0}$, $\widehat{\boldsymbol{\beta}}_{T}$
converges $\boldsymbol{\beta}_{T}^{0}$, and $\widehat{R}_{T}(\cdot)$ converges
to $R_{T}^{0}(\cdot)$ uniformly on $[0,\xi]$.

###  C.2 Proof of Theorem 2

By Theorem 3.3.1 of van der Vaart and Wellner, (1996),
$\sqrt{n}(\widehat{\boldsymbol{\theta}}_{1}-\boldsymbol{\theta}_{1}^{0})$ can
be expressed as

$\sqrt{n}(\widehat{\boldsymbol{\theta}}_{1}-\boldsymbol{\theta}_{1}^{0})=\mathcal{I}_{1,n}^{-1}(\boldsymbol{\theta}_{1}^{0},\widehat{\boldsymbol{\theta}}_{D})\left\\{n^{-1/2}\sum_{i=1}^{n}\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},\widehat{\boldsymbol{\theta}}_{D})\right\\}+o_{p}(1),$

where

$\mathcal{I}_{1,n}(\boldsymbol{\theta}_{1}^{0},\widehat{\boldsymbol{\theta}}_{D})=-n^{-1}\sum_{i=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}_{1}\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},\widehat{\boldsymbol{\theta}}_{D})$
(10)

with
$\ddot{\ell}_{\boldsymbol{\theta}_{1}\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},\widehat{\boldsymbol{\theta}}_{D})=\partial^{2}\ell_{i}(\boldsymbol{\theta})/\partial\boldsymbol{\theta}_{1}\partial\boldsymbol{\theta}_{1}^{{\sf\scriptscriptstyle{T}}}$.
By C6 and the consistency of $\widehat{\boldsymbol{\theta}}_{D}$,

$\left\|\mathcal{I}_{1,n}(\boldsymbol{\theta}_{1}^{0},\widehat{\boldsymbol{\theta}}_{D})-\mathbb{E}\left\\{-\ddot{\ell}_{\boldsymbol{\theta}_{1}\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}^{0})\right\\}\right\|=o_{p}(1).$

Thus,

$\sqrt{n}(\widehat{\boldsymbol{\theta}}_{1}-\boldsymbol{\theta}_{1}^{0})=\left(\mathcal{I}_{1}^{0}\right)^{-1}\left\\{n^{-1/2}\sum_{i=1}^{n}\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},\widehat{\boldsymbol{\theta}}_{D})\right\\}+o_{p}(1),$

where
$\mathcal{I}_{1}^{0}=\mathbb{E}\\{-\ddot{\ell}_{\boldsymbol{\theta}_{1}\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}^{0})\\}$.
Since
$\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},\widehat{\boldsymbol{\theta}}_{D})$
is a function of $\widehat{\boldsymbol{\theta}}_{D}$ only through
$\widehat{U}_{D,i}$, we can write
$\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},\widehat{\boldsymbol{\theta}}_{D})=\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},\widehat{U}_{D,i})$.
By the mean value theorem, we have

$\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},\widehat{U}_{D,i})=\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},U_{D,i}^{0})+\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},i}(\boldsymbol{\theta}_{1}^{0},U_{D,i}^{0})(\widehat{U}_{D,i}-U_{D,i}^{0})+o_{p}(1),$

where
$\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},i}(\boldsymbol{\theta}_{1}^{0},U_{D,i}^{0})=\partial\Psi_{1,i}(\boldsymbol{\theta}_{1}^{0})/\partial
u_{D}=\partial^{2}\ell_{i}(\boldsymbol{\theta}^{0})/\partial
u_{D}\partial\boldsymbol{\theta}_{1}$. Thus, we can express

$\displaystyle\sqrt{n}(\widehat{\boldsymbol{\theta}}_{1}-\boldsymbol{\theta}_{1}^{0})=\left(\mathcal{I}_{1}^{0}\right)^{-1}n^{-1/2}\sum_{i=1}^{n}\left\\{\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}_{1}^{0},U_{D,i}^{0})+\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},i}(\boldsymbol{\theta}_{1}^{0},U_{D,i}^{0})(\widehat{U}_{D,i}-U_{D,i}^{0})\right\\}+o_{p}(1).$

By the expansion of $\sqrt{n}(\widehat{U}_{D,i}-U_{D,i}^{0})$ in (9), we can
express
$\sqrt{n}(\widehat{\boldsymbol{\theta}}_{1}-\boldsymbol{\theta}_{1}^{0})=n^{-1/2}\sum_{i=1}^{n}\psi_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}^{0})\allowbreak+o_{p}(1)$,
where

$\psi_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}^{0})=\left(\mathcal{I}_{1}^{0}\right)^{-1}\left\\{\dot{\ell}_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}^{0})+n^{-1}\sum_{k=1}^{n}\ddot{\ell}_{\boldsymbol{\theta}_{1}u_{D},k}(\boldsymbol{\theta}^{0})\psi_{u_{D,k},i}(\boldsymbol{\theta}_{D}^{0})\right\\}.$

Thus,
$\sqrt{n}\\{F_{1}(\widehat{\boldsymbol{\theta}}_{1})-F_{1}(\boldsymbol{\theta}_{1}^{0})\\}$
converges to a zero-mean Gaussian process with variance-covariance matrix
$\dot{F}_{1}(\boldsymbol{\theta}_{1}^{0})^{{\sf\scriptscriptstyle{T}}}\boldsymbol{\Sigma}_{1}^{0}\dot{F}_{1}(\boldsymbol{\theta}_{1}^{0})$,
where
$\boldsymbol{\Sigma}_{1}^{0}=\mathrm{Var}\\{\psi_{\boldsymbol{\theta}_{1},i}(\boldsymbol{\theta}^{0})\\}$.

| $n=200$ | $n=400$
---|---|---
$\tau$ | $\boldsymbol{\beta}_{T}$ | Method | BIAS | ESD | rMSE | BIAS | ESD | rMSE
| $\widehat{\beta}_{T,1}$
0.4 | (1,1) | MLE | 0.005 | 0.160 | 0.160 | 0.001 | 0.113 | 0.113
PMLE | 0.006 | 0.162 | 0.162 | 0.002 | 0.115 | 0.115
(1,0.5) | MLE | 0.006 | 0.163 | 0.163 | 0.001 | 0.114 | 0.114
PMLE | 0.007 | 0.165 | 0.166 | 0.002 | 0.116 | 0.116
0.6 | (1,1) | MLE | 0.001 | 0.156 | 0.156 | 0.002 | 0.111 | 0.111
PMLE | 0.002 | 0.163 | 0.163 | 0.001 | 0.117 | 0.117
(1,0.5) | MLE | 0.000 | 0.158 | 0.158 | 0.001 | 0.113 | 0.113
PMLE | 0.003 | 0.164 | 0.164 | 0.002 | 0.118 | 0.118
0.8 | (1,1) | MLE | 0.007 | 0.156 | 0.156 | 0.006 | 0.111 | 0.112
PMLE | 0.005 | 0.168 | 0.168 | 0.002 | 0.121 | 0.121
(1,0.5) | MLE | 0.005 | 0.157 | 0.157 | 0.004 | 0.112 | 0.113
PMLE | 0.006 | 0.169 | 0.169 | 0.003 | 0.121 | 0.121
| $\widehat{\beta}_{T,2}$
0.4 | (1,1) | MLE | 0.012 | 0.203 | 0.203 | 0.003 | 0.148 | 0.148
PMLE | 0.009 | 0.204 | 0.205 | 0.002 | 0.150 | 0.150
(1,0.5) | MLE | 0.024 | 0.393 | 0.394 | 0.009 | 0.285 | 0.285
PMLE | 0.020 | 0.396 | 0.397 | 0.007 | 0.289 | 0.289
0.6 | (1,1) | MLE | 0.016 | 0.198 | 0.198 | 0.008 | 0.142 | 0.142
PMLE | 0.009 | 0.204 | 0.204 | 0.005 | 0.148 | 0.148
(1,0.5) | MLE | 0.025 | 0.383 | 0.383 | 0.017 | 0.274 | 0.274
PMLE | 0.015 | 0.394 | 0.394 | 0.012 | 0.285 | 0.285
0.8 | (1,1) | MLE | 0.024 | 0.196 | 0.197 | 0.011 | 0.139 | 0.139
PMLE | 0.006 | 0.211 | 0.211 | 0.003 | 0.151 | 0.151
(1,0.5) | MLE | 0.033 | 0.382 | 0.383 | 0.018 | 0.269 | 0.270
PMLE | 0.011 | 0.409 | 0.409 | 0.009 | 0.291 | 0.291
| $\widehat{\alpha}$
0.4 | (1,1) | MLE | 0.036 | 0.075 | 0.083 | 0.018 | 0.051 | 0.054
PMLE | 0.023 | 0.073 | 0.077 | 0.011 | 0.051 | 0.052
(1,0.5) | MLE | 0.037 | 0.078 | 0.086 | 0.019 | 0.053 | 0.056
PMLE | 0.024 | 0.075 | 0.079 | 0.012 | 0.052 | 0.054
0.6 | (1,1) | MLE | 0.054 | 0.087 | 0.103 | 0.027 | 0.057 | 0.064
PMLE | 0.026 | 0.083 | 0.087 | 0.013 | 0.057 | 0.059
(1,0.5) | MLE | 0.055 | 0.089 | 0.105 | 0.028 | 0.059 | 0.065
PMLE | 0.027 | 0.084 | 0.088 | 0.013 | 0.058 | 0.060
0.8 | (1,1) | MLE | 0.096 | 0.105 | 0.142 | 0.049 | 0.065 | 0.081
PMLE | 0.020 | 0.093 | 0.095 | 0.010 | 0.062 | 0.063
(1,0.5) | MLE | 0.094 | 0.106 | 0.142 | 0.049 | 0.065 | 0.082
PMLE | 0.020 | 0.093 | 0.096 | 0.010 | 0.063 | 0.064
Table 1: Simulation results for the estimates under the Gumbel copula. BIAS: relative bias; ESD: relative empirical standard deviation; rMSE: relative root mean square error. | $n=200$ | $n=400$
---|---|---
$\tau$ | $\boldsymbol{\beta}_{T}$ | Method | BIAS | ESD | rMSE | BIAS | ESD | rMSE
| $\widehat{\beta}_{T,1}$
0.4 | (1,1) | MLE | 0.014 | 0.154 | 0.155 | 0.008 | 0.108 | 0.108
PMLE | 0.014 | 0.160 | 0.161 | 0.008 | 0.112 | 0.112
(1,0.5) | MLE | 0.016 | 0.160 | 0.161 | 0.008 | 0.112 | 0.112
PMLE | 0.016 | 0.166 | 0.167 | 0.008 | 0.115 | 0.115
0.6 | (1,1) | MLE | 0.018 | 0.143 | 0.144 | 0.011 | 0.101 | 0.101
PMLE | 0.018 | 0.161 | 0.162 | 0.011 | 0.114 | 0.115
(1,0.5) | MLE | 0.022 | 0.147 | 0.149 | 0.011 | 0.103 | 0.103
PMLE | 0.021 | 0.165 | 0.167 | 0.013 | 0.116 | 0.117
0.8 | (1,1) | MLE | 0.034 | 0.133 | 0.137 | 0.016 | 0.091 | 0.093
PMLE | 0.027 | 0.173 | 0.175 | 0.016 | 0.122 | 0.123
(1,0.5) | MLE | 0.035 | 0.135 | 0.140 | 0.016 | 0.091 | 0.093
PMLE | 0.029 | 0.174 | 0.177 | 0.018 | 0.123 | 0.124
| $\widehat{\beta}_{T,2}$
0.4 | (1,1) | MLE | 0.009 | 0.197 | 0.197 | 0.002 | 0.145 | 0.145
PMLE | 0.003 | 0.207 | 0.207 | 0.008 | 0.150 | 0.150
(1,0.5) | MLE | 0.026 | 0.372 | 0.373 | 0.007 | 0.275 | 0.276
PMLE | 0.008 | 0.392 | 0.393 | 0.002 | 0.284 | 0.284
0.6 | (1,1) | MLE | 0.008 | 0.177 | 0.177 | 0.001 | 0.125 | 0.125
PMLE | 0.012 | 0.210 | 0.210 | 0.012 | 0.147 | 0.148
(1,0.5) | MLE | 0.023 | 0.332 | 0.333 | 0.012 | 0.236 | 0.236
PMLE | 0.005 | 0.392 | 0.392 | 0.007 | 0.279 | 0.279
0.8 | (1,1) | MLE | 0.001 | 0.161 | 0.161 | 0.002 | 0.111 | 0.111
PMLE | 0.025 | 0.223 | 0.225 | 0.018 | 0.156 | 0.157
(1,0.5) | MLE | 0.011 | 0.303 | 0.304 | 0.012 | 0.208 | 0.209
PMLE | 0.019 | 0.413 | 0.414 | 0.014 | 0.291 | 0.292
| $\widehat{\alpha}$
0.4 | (1,1) | MLE | 0.063 | 0.186 | 0.196 | 0.025 | 0.124 | 0.127
PMLE | 0.026 | 0.184 | 0.186 | 0.007 | 0.125 | 0.125
(1,0.5) | MLE | 0.061 | 0.191 | 0.201 | 0.026 | 0.127 | 0.129
PMLE | 0.025 | 0.189 | 0.190 | 0.007 | 0.127 | 0.127
0.6 | (1,1) | MLE | 0.067 | 0.148 | 0.162 | 0.028 | 0.094 | 0.099
PMLE | 0.011 | 0.146 | 0.146 | 0.001 | 0.098 | 0.098
(1,0.5) | MLE | 0.065 | 0.151 | 0.165 | 0.027 | 0.095 | 0.099
PMLE | 0.010 | 0.149 | 0.150 | 0.001 | 0.098 | 0.098
0.8 | (1,1) | MLE | 0.101 | 0.138 | 0.171 | 0.046 | 0.080 | 0.092
PMLE | 0.011 | 0.125 | 0.126 | 0.011 | 0.080 | 0.081
(1,0.5) | MLE | 0.098 | 0.139 | 0.170 | 0.045 | 0.081 | 0.093
PMLE | 0.013 | 0.127 | 0.127 | 0.012 | 0.081 | 0.081
Table 2: Simulation results for the estimates under the Clayton copula. BIAS:
relative bias; ESD: relative empirical standard deviation; rMSE: relative root
mean square error.

| $\boldsymbol{\beta}_{T}=(1,0.5)$ | $\boldsymbol{\beta}_{T}=(1,1)$
---|---|---
| $\widehat{\beta}_{T,1}$ | $\widehat{\beta}_{T,2}$ | $\widehat{\alpha}$ | $\widehat{\beta}_{T,1}$ | $\widehat{\beta}_{T,2}$ | $\widehat{\alpha}$
$\tau$ | $n$ | MLE | PMLE | MLE | PMLE | MLE | PMLE | MLE | PMLE | MLE | PMLE | MLE | PMLE
Gumbel
0.4 | 200 | 1.002 | 0.997 | 1.041 | 1.030 | 0.974 | 0.979 | 1.000 | 0.994 | 1.028 | 1.018 | 0.971 | 0.978
400 | 0.993 | 0.996 | 1.069 | 1.069 | 0.978 | 0.988 | 1.001 | 1.002 | 1.065 | 1.065 | 0.977 | 0.984
0.6 | 200 | 1.000 | 0.988 | 1.045 | 1.027 | 0.997 | 0.990 | 0.998 | 0.985 | 1.034 | 1.016 | 0.990 | 0.988
400 | 1.007 | 1.010 | 1.062 | 1.063 | 0.980 | 0.992 | 1.007 | 1.010 | 1.055 | 1.056 | 0.976 | 0.988
0.8 | 200 | 1.001 | 0.988 | 1.054 | 1.040 | 1.013 | 0.987 | 0.999 | 0.985 | 1.039 | 1.026 | 1.004 | 0.981
400 | 1.010 | 1.006 | 1.049 | 1.057 | 0.977 | 0.978 | 1.009 | 1.007 | 1.043 | 1.051 | 0.972 | 0.970
Clayton
0.4 | 200 | 1.018 | 1.004 | 1.038 | 1.031 | 1.005 | 1.012 | 1.015 | 1.001 | 1.034 | 1.027 | 0.999 | 1.006
400 | 1.013 | 0.993 | 1.099 | 1.074 | 0.980 | 0.982 | 1.014 | 0.997 | 1.085 | 1.067 | 0.983 | 0.985
0.6 | 200 | 1.033 | 0.997 | 1.048 | 1.042 | 1.020 | 1.019 | 1.030 | 0.992 | 1.037 | 1.045 | 1.010 | 1.008
400 | 1.033 | 0.993 | 1.073 | 1.066 | 0.962 | 0.974 | 1.042 | 0.999 | 1.060 | 1.055 | 0.970 | 0.980
0.8 | 200 | 1.005 | 0.997 | 1.027 | 1.051 | 1.017 | 0.987 | 0.996 | 0.996 | 1.024 | 1.062 | 1.017 | 0.977
400 | 1.020 | 0.993 | 1.048 | 1.059 | 0.938 | 0.927 | 1.027 | 0.997 | 1.042 | 1.059 | 0.933 | 0.923

Table 3: Simulation results for the ratio of the ESD to ASE. ESD: empirical standard deviation; ASE: average standard error. Parameter | Method | BIAS | ESD | MSE
---|---|---|---|---
True: Gumbel | GC | GG | GC | GG | GC | GG
$\beta_{D,1}$ | MLE | 0.070 | 0.003 | 0.108 | 0.111 | 0.129 | 0.111
PMLE | 0.002 | 0.002 | 0.118 | 0.118 | 0.118 | 0.118
$\beta_{D,2}$ | MLE | 0.019 | 0.002 | 0.122 | 0.134 | 0.124 | 0.134
PMLE | 0.001 | 0.001 | 0.144 | 0.144 | 0.144 | 0.144
$\beta_{T,1}$ | MLE | 0.195 | 0.001 | 0.116 | 0.113 | 0.226 | 0.113
PMLE | 0.078 | 0.002 | 0.115 | 0.118 | 0.139 | 0.118
$\beta_{T,2}$ | MLE | 0.104 | 0.008 | 0.130 | 0.137 | 0.166 | 0.137
PMLE | 0.045 | 0.006 | 0.137 | 0.143 | 0.144 | 0.143
True: Clayton | CG | CC | CG | CC | CG | CC
$\beta_{D,1}$ | MLE | 0.053 | 0.005 | 0.123 | 0.095 | 0.133 | 0.095
PMLE | 0.003 | 0.003 | 0.117 | 0.117 | 0.117 | 0.117
$\beta_{D,2}$ | MLE | 0.021 | 0.013 | 0.149 | 0.112 | 0.151 | 0.113
PMLE | 0.001 | 0.001 | 0.145 | 0.145 | 0.145 | 0.145
$\beta_{T,1}$ | MLE | 0.041 | 0.011 | 0.123 | 0.103 | 0.129 | 0.103
PMLE | 0.041 | 0.013 | 0.125 | 0.116 | 0.131 | 0.117
$\beta_{T,2}$ | MLE | 0.013 | 0.006 | 0.152 | 0.118 | 0.152 | 0.118
PMLE | 0.018 | 0.003 | 0.154 | 0.140 | 0.155 | 0.140
Table 4: Simulation results for copula misspecification. GC: true Gumbel - fit Clayton; GG: Gumbel - Gumbel; GC: Gumbel - Clayton; CG: Clayton - Gumbel; CC: Clayton - Clayton; CG: Clayton - Gumbel. BIAS: bias; ESD: empirical standard deviation; rMSE: root mean square error. | $n=200$ | $n=400$
---|---|---
$\tau$ | $\alpha^{*}$ | $\boldsymbol{\beta}_{T}$ | Method | BIAS | ESD | rMSE | BIAS | ESD | rMSE
0.4 | 0.6 | (1,1) | PMLE | 0.011 | 0.042 | 0.043 | 0.005 | 0.030 | 0.030
Zhu et al., (2021) | 0.034 | 0.059 | 0.068 | 0.028 | 0.058 | 0.064
(1,0.5) | PMLE | 0.011 | 0.043 | 0.044 | 0.006 | 0.031 | 0.031
Zhu et al., (2021) | 0.032 | 0.061 | 0.069 | 0.029 | 0.059 | 0.066
0.6 | 0.4 | (1,1) | PMLE | 0.008 | 0.031 | 0.032 | 0.004 | 0.022 | 0.023
Zhu et al., (2021) | 0.058 | 0.064 | 0.086 | 0.041 | 0.057 | 0.070
(1,0.5) | PMLE | 0.008 | 0.032 | 0.033 | 0.004 | 0.023 | 0.023
Zhu et al., (2021) | 0.051 | 0.060 | 0.079 | 0.039 | 0.056 | 0.068
0.8 | 0.2 | (1,1) | PMLE | 0.002 | 0.018 | 0.018 | 0.001 | 0.012 | 0.012
Zhu et al., (2021) | 0.078 | 0.056 | 0.096 | 0.057 | 0.042 | 0.071
(1,0.5) | PMLE | 0.002 | 0.018 | 0.018 | 0.001 | 0.012 | 0.012
Zhu et al., (2021) | 0.068 | 0.048 | 0.083 | 0.055 | 0.042 | 0.069
Table 5: Simulation results comparing the proposed PMLE and the method in Zhu
et al., (2021). BIAS: bias; ESD: empirical standard deviation; rMSE: root mean
square error. $\alpha^{*}$ is the copula parameter of Gumbel copula given by
$\mathbbm{C}(u_{1},u_{2};\alpha^{*})=\exp\left(-\left[\\{-\log(u_{1})\\}^{1/\alpha^{*}}+\\{-\log(u_{2})\\}^{1/\alpha^{*}}\right]^{\alpha^{*}}\right)$.

Covariates | Gumbel | Clayton | Frank
---|---|---|---
| PMLE | MLE | PMLE | MLE | PMLE | MLE
Regression coefficients for leukemia relapse
AML high | 1.239 (0.317) | 1.147 (0.313) | 1.168 (0.311) | 1.116 (0.306) | 1.137 (0.302) | 1.032 (0.305)
ALL | 0.854 (0.345) | 0.764 (0.325) | 0.710 (0.324) | 0.669 (0.320) | 0.716 (0.318) | 0.645 (0.318)
Regression coefficients for death
AML high | 1.022 (0.276) | 0.953 (0.281) | 1.022 (0.276) | 0.977 (0.271) | 1.022 (0.276) | 0.905 (0.272)
ALL | 0.611 (0.285) | 0.553 (0.284) | 0.611 (0.285) | 0.577 (0.280) | 0.611 (0.285) | 0.535 (0.279)
Kendall’s $\tau$ between event times
AML low | 0.679 (0.140) | 0.709 (0.122) | 0.814 (0.084) | 0.821 (0.079) | 0.737 (0.105) | 0.738 (0.106)
AML high | 0.726 (0.090) | 0.769 (0.089) | 0.769 (0.108) | 0.777 (0.102) | 0.760 (0.094) | 0.770 (0.095)
ALL | 0.686 (0.090) | 0.726 (0.079) | 0.767 (0.100) | 0.773 (0.093) | 0.721 (0.094) | 0.730 (0.095)
Table 6: Estimation results on the BMT data using the PMLE and MLE methods.
Standard errors are given in parentheses.

Figure 1: PMLE and MLE of relapse-free survival probabilities
$\mathbb{P}(T>t)$ for each disease group.