A Simulation Study on Comparing General Class of Semiparametric Transformation Models for Survival Outcome with Time-Varying Coefficients and Covariates

Yemane Hailu Fissuh^{1,2}^{*},
Tsegay Giday Woldu^{3},
Idriss Abdelmajid Idriss Ahmed^{1},
Abebe Zewdie Kebebe^{4}

Show more

1. Introduction

In many experimental and observational studies such as randomized clinical trials, agricultural experiments, and engineering and industrial production commonly we obtain time-to-end outcomes so-called survival time or failure time. In biomedical researches, the main concern is usually on the survival time, which is a time from defined origin until the defined endpoint or outcome [1] . The survival data have missing value raised through the censoring mechanisms. Censoring is the problem of not finding the exact time of an event during the experimental or observational studies, which makes the analysis much more complex.

Central to the entire discipline of survival analysis, mostly right censoring exists. Besides, a time-varying covariate is a classical problem in modeling survival time. The semiparametric transformation models which have been attracted by several authors have been an important concept in the study of right censored survival time. The another important concept in analysing survival data is proportionality assumption. Sometimes, in our experimental study, we have no warrantee of the fulfillment of this assumption. Because the effect of covariate may vary over time breaking the proportionality assumption for Cox proportional hazards model of [2] . In this situation, we need to consider the time-varying coefficient to our model. Due to this, the time-dependent effect and time-dependent covariates have been given attentions these days. Generally, someone may need to extend this model to more general model that can incorporate both time-varying covariate and time-varying effect. Thus the combination brings more general version.

A key role of semiparametric transformation models (STM) is that the model provides a framework for deriving the effect of time-varying covariates and the effect of time-varying coefficients on failure time. In this model, since the model consists of different special cases inside, the failure of proportionality assumption might not be much problem.

The remaining part of this paper is organized as follows. Section 2 introduces the methods and model framework which are going to be used in the whole paper and proposes a modified estimating equation for robust semiparametric transformation models. Section 3 presents a large sample theory and regularity conditions for the consistency and asymptotic properties of the proposed estimators. Section 4 devotes simulation studies to check the performance of the proposed techniques. Finally, the conclusion is presented in Section 5.

2. Methods and Model Framework

Here we start with some basic notations that are used throughout this paper. $T\ge 0$ which is $T\in {\mathbb{R}}^{+}$ be the survival time of interest, and $X\left(t\right)$ be $p\times 1$ vectors of possibly time-varying covariates, where p is number covariates included in the model.

Where the covariate is allowed to vary over time, possibly the furthermost instant tactic is to use the step-function as follows.

$X\left(t\right)=\{\begin{array}{cc}1,& \text{if}\text{\hspace{0.17em}}t>{T}_{r},\\ 0,& t\le {T}_{r},\end{array}$ (1)

where ${T}_{r}$ is transaction time for change.

Whenever the covariate only changes once at fixed time point and do not change after that, the step function is used. However, in some situations it is common to have covariate that change over time continuously and frequently at a time with the only requirement that the intervals of the observation need not be contiguous. Therefore, in this situation a simple way to code time-dependent covariates is using intervals of time and recorded in to two columns as the start, stop or time 1, time 2, entry, end and so forth. The “tmerge” package in R can do this arrangement in the survival library.

When the censoring time is denoted by C, the failure or censor time represented by $\stackrel{\u02dc}{Y}$ is the minimum of failure time of censoring time and failure time; i.e.; $\stackrel{\u02dc}{Y}=\mathrm{min}\left(T,C\right)$ . We write $\delta =\mathbb{I}\left(T\le C\right)$ for the event indicator. Finally, the summarized n independent random vectors of observations are formulated as

${O}_{i}=\left\{{\stackrel{\u02dc}{Y}}_{i}\mathrm{,}{\delta}_{i}\mathrm{,}{X}_{i}\left(t\right)\right\}\mathrm{.}$ (2)

2.1. The Semiparametric Transformation Models

The flexibility extended general class of semiparametric transformation models with the effect of time-varying coefficients is formulated

$\mathcal{H}\left(t\right)=-{\vartheta}^{\prime}\left(t\right)X+\epsilon \mathrm{,}\epsilon =\mathrm{log}\u03f5\mathrm{,}$ (3)

where X is a set of covariates, the set of time-varying regression coefficients or parameters $\vartheta \left(t\right)=\mathrm{log}\left(T\right)$ , where $\mathrm{log}\left(T\right)$ is a natural logarithm or logarithm of base $\text{e}=2.71828\cdots $ and the unspecified continuously differentiable monotone arbitrary transformation function $\mathcal{H}\left(\mathrm{.}\right)=\mathrm{log}{\Lambda}_{0}\left(\mathrm{.}\right)$ satisfying $\mathcal{H}\left(0\right)=-\infty $ , are unknown and the extraneous random error term $\epsilon =\mathrm{log}\u03f5\perp |X$ comes from unrestricted well-known parametric distributions ${F}_{\epsilon}$ ; see [3] . The unspecified monotonically increasing transformation function ${\Lambda}_{0}\left(t\right)={\displaystyle {\int}_{0}^{t}{\lambda}_{0}\left(u\right)\text{d}u}$ is cumulative baseline hazard function satisfying ${\Lambda}_{0}\left(0\right)=0$ . Then the conditional distribution is formulated as

${F}_{T\mathrm{|}X}\left(t\mathrm{|}x\right)=Pr\left(T\le t\mathrm{|}X=x\right)={F}_{\epsilon}\left(\mathcal{H}\left(t\right)+{\vartheta}^{\prime}\left(t\right)X\right)\mathrm{.}$ (4)

However, the model (3) does not applicable for time-varying covariate. Then, with the extension of time-varying covariates, the special cases of the transformation models consider proportional hazards (PH) model and proportional odds (PO) model. These special models are based on the given distribution of random error term $\epsilon $ corresponds to extreme value distribution and the standard logistic distribution respectively [4] [5] [6] .

Let ${\mathcal{N}}_{i}\left(t\right)$ be the counting process recording the number of events that have occurred by time t and let $X\left(t\right)$ be a set of predictors which contains a vector of possibly time-varying covariates. We specify that the cumulative intensity function for ${\mathcal{N}}_{i}\left(t\right)$ conditional on $X\left(t\right)\mathrm{:}X\left(u\right)\mathrm{;}u\le t\}$ and therefore, equivalent formulation of model (3) can be expressed as

$\Lambda \left\{t\mathrm{|}X\left(t\right)=x\left(t\right)\right\}=\Phi \left({\displaystyle {\int}_{0}^{t}}\text{\hspace{0.05em}}\text{d}{\Lambda}_{\epsilon}\left\{\mathcal{H}\left(u\right)+{\vartheta}^{\prime}\left(u\right)X\left(u\right)\right\}\right)\mathrm{,}$ (5)

where $\Lambda \left\{t|X\left(t\right)=x\left(t\right)\right\}$ is subject-specific cumulative hazards function under completely specified continuous monotonically increasing transformation function $\Phi \left(\mathrm{.}\right)$ satisfying for $\Phi \left(0\right)=0$ and $\Phi \left(\infty \right)=\infty $ . Here the independent identically distributed (i.i.d) random variable $\u03f5$ with known distribution is unobservable positive noise associated to random biological features.

For strictly increasing transformation function $\Phi \left(\mathrm{.}\right)$ , the class of Box-Cox transformations which was recently used by [7] is also considered here. For the two special cases of transformation model classes namely Proportional hazards (PH) and proportional odds (PO) models, we reflect on the Box-Cox transformation functions

$\Phi \left(x\right)=\{\begin{array}{l}\frac{{\left(1+x\right)}^{\rho}-1}{\rho},\mathrm{}\rho >0,\hfill \\ \mathrm{log}\left(1+x\right),\mathrm{}\rho =0,\hfill \end{array}$ (6)

and the class of logarithmic transformation

$\Phi \left(x\right)=\{\begin{array}{l}\frac{\mathrm{log}\left(1+rx\right)}{r},\mathrm{}r>0,\hfill \\ x,\mathrm{}\text{\hspace{0.05em}}\text{\hspace{0.05em}}r=0.\hfill \end{array}$ (7)

Therefore, the choice of $\Phi \left(x\right)=x$ when either $\rho =1$ or $r=0$ , the special case of transformation model indeed yields PH model for survival data. Equivalently, the choice of $\Phi \left(x\right)=\mathrm{log}\left(1+x\right)$ for either $\rho =0$ or $r=1$ , the special case of transformation model indeed yields PO model for survival data.

Remark: Specifying the function $\Phi $ while leaving the function ${\Lambda}_{0}$ unspecified is equivalent to specifying the distribution of $\epsilon $ while leaving the function $\mathcal{H}$ unspecified. Non-identifiability arises if both $\Phi $ and ${\Lambda}_{0}$ (or both $\mathcal{H}$ and $\epsilon $ ) are unspecified and $\vartheta =0$ ( [3] , p. 169) which was quoted by [8] .

The Modified Estimating Equations

Before developing estimating equations, let us impose on the following two unignorable assumptions.

Assumption 1: The parameter space of $\Theta $ of $\vartheta $ is bounded open subset of ${\mathbb{R}}^{k}$ .

Assumption 2: The random variable censoring time C is the independent of random variable time of failure or event T given possibly time varying covariates $X\left(t\right)$ i.e.; $C\perp \left|T\right|X\left(t\right)$ .

To develop the estimating equations to estimate the unknown parameter $\vartheta $ and unknown strictly increasing monotonic function $\mathcal{H}\left(u\right)$ , estimating equations of [5] which has been lately used by several authors for example [9] [10] [11] and [12] is modified for the effect of time-varying coeffcients and time-varying covariates.

In this paper we suppose, $0\le {t}_{1}\le {t}_{2}\le \cdots \le {t}_{r}\le \infty $ for the r failure times among the n observations. Furthermore, we suppose ${\lambda}_{\epsilon}(.)$ and ${\Lambda}_{\epsilon}(.)$ be the known hazard and cumulative hazard functions of $\epsilon $ , respectively. Let us propose the true values of $\mathcal{H}\left(t\right)$ and $\vartheta \left(t\right)$ denoted by ${\mathcal{H}}_{0}(.)$ and ${\mathcal{H}}_{0}(.)$ respectively. Therefore, following the usual counting process notation, let

$Y\left(t\right)=I\left({\stackrel{\u02dc}{Y}}_{i}\ge t\right),{N}_{i}\left(t\right)={\delta}_{i}I\left({\stackrel{\u02dc}{Y}}_{i}\le t,\delta =1\right),$ (8)

are an at-risk indicator process and the distinct ordered uncensored failure times ${y}_{1}\mathrm{,}\cdots \mathrm{,}{y}_{r}$ respectively. Suppose $\left\{{Y}_{i}\left(t\right)\mathrm{,}{N}_{i}\left(t\right)\mathrm{,}{M}_{i}\left(t\right)\right\}$ be the sample analogues of $\left\{Y\left(t\right)\mathrm{,}N\left(t\right)\mathrm{,}M\left(t\right)\right\}$ . Thus, the martingale decomposition can minimize the complexity of the estimation of equations by constructing the following easily tractable formula.

${M}_{i}\left(t\right)={N}_{i}\left(t\right)-\Phi \left({\displaystyle {\int}_{0}^{t}}\text{\hspace{0.05em}}{Y}_{i}\left(u\right)\text{d}{\Lambda}_{\epsilon}\left\{{\mathcal{H}}_{0}\left(u\right)+{{\vartheta}^{\prime}}_{0}\left(u\right)X\left(u\right)\right\}\right)\mathrm{,}$ (9)

for complete σ-field $\mathcal{F}$ since

$\begin{array}{l}E\left[\text{d}{N}_{i}\left(t\mathrm{|}\mathcal{F}\left(t-\right)\right)\right]\\ ={Y}_{i}\left(t\right){\lambda}_{\epsilon}\left\{{\stackrel{^}{\vartheta}}^{\prime}\left(t\right){X}_{i}\left(t\right)+\stackrel{^}{\mathcal{H}}\left(t\right)\right\}\stackrel{\dot{}}{\Phi}\left({\displaystyle {\int}_{0}^{t}}{Y}_{i}\left(u\right)\text{d}{\Lambda}_{\epsilon}\left\{{\mathcal{H}}_{0}\left(u\right)+{{\vartheta}^{\prime}}_{0}\left(u\right){X}_{i}\left(u\right)\right\}\right)\mathrm{,}\end{array}$ (10)

where $\stackrel{\dot{}}{\Phi}\left(\mathrm{.}\right)=\frac{\text{d}}{\text{d}t}\Phi \left(\mathrm{.}\right)$ , $\stackrel{^}{\vartheta}\left(t\right)$ and $\stackrel{^}{\mathcal{H}}\left(t\right)$ denote the estimators of $\vartheta \left(t\right)$ and $\mathcal{H}\left(t\right)$ and the mean of a martingale process with respect to $\mathcal{F}$ is zero.

Lemma 1: The mean of the derivative of regular martingale process is zero.

$E\left(\text{d}{M}_{i}\left(t\right)\right)=0.$ (11)

Thus, slightly modified estimating equations of [5] are proposed by making possibly time-varying covariate under consideration. The two modified estimating equations are

${Q}_{1}\left(\vartheta \left(t\right)\mathrm{,}\mathcal{H}\left(\mathrm{.}\right)\right)={\displaystyle \underset{i=1}{\overset{n}{\sum}}}{\displaystyle {\int}_{0}^{t}}{X}_{i}\left(u\right)\left[\text{d}{N}_{i}\left(u\right)-{\lambda}_{T\mathrm{|}{X}_{i}\left(t\right)}\left\{{\vartheta}^{\prime}\left(u\right){X}_{i}\left(u\right)+\mathcal{H}\left(u\right)\right\}\right]=\mathrm{0,}$ (12)

${Q}_{2}\left(\vartheta \left(t\right)\mathrm{,}\mathcal{H}\left(\mathrm{.}\right)\right)={\displaystyle \underset{i=1}{\overset{n}{\sum}}}\left[\text{d}{N}_{i}\left(t\right)-{\lambda}_{T\mathrm{|}X}{{}_{{}_{i}}}_{\left(t\right)}\left\{{\vartheta}^{\prime}\left(t\right){X}_{i}\left(t\right)+\mathcal{H}\left(t\right)\right\}\right]=\mathrm{0,}\left(t\ge 0\right)\mathrm{,}$ (13)

where

$\begin{array}{l}{\lambda}_{T\mathrm{|}X}{{}_{{}_{i}}}_{\left(t\right)}\left\{{\vartheta}^{\prime}\left(t\right){X}_{i}\left(t\right)+\mathcal{H}\left(t\right)\right\}\\ ={Y}_{i}\left(t\right){\lambda}_{\epsilon}\left\{{\vartheta}^{\prime}\left(u\right){X}_{i}\left(u\right)+\mathcal{H}\left(u\right)\right\}\stackrel{\dot{}}{\Phi}\left({\displaystyle {\int}_{0}^{t}}{Y}_{i}\left(u\right)\text{d}{\Lambda}_{\epsilon}\left\{{\vartheta}^{\prime}\left(u\right){X}_{i}\left(u\right)+\mathcal{H}\left(u\right)\right\}\right)\mathrm{,}\end{array}$ (14)

is the intensity function for ${N}_{i}\left(t\right)$ and $\mathcal{H}\left(\mathrm{.}\right)$ is nondecreasing function satisfying $\mathcal{H}\left(0\right)=-\infty $ . Therefore, this requirement in turn ensures that for any finite number k, ${\Lambda}_{\epsilon}\left\{k+\mathcal{H}\left(0\right)\right\}=0$ .

For the special case when we assume the Cox’s proportional hazards model of [5] in which $\lambda \left(t\right)=\Lambda \left(t\right)=\mathrm{exp}\left(t\right)$ , while $\Phi \left(x\right)=x$ , and

$\text{d}\left(\mathrm{exp}\left(\mathcal{H}\left(t\right)\right)\right)+{\vartheta}^{\prime}\left(t\right)\ast \mathrm{exp}\left(\mathcal{H}\left(t\right)\right)=\frac{{\displaystyle {\sum}_{i=1}^{n}}\text{\hspace{0.05em}}\text{d}{N}_{i}\left(t\right)}{{\displaystyle {\sum}_{i=1}^{n}}{Y}_{i}\left(t\right)\ast \mathrm{exp}\left({\vartheta}^{\prime}\left(t\right){X}_{i}\left(t\right)\right)}$ ,

therefore, by plugging this in (12) we simply obtain

$\underset{i=1}{\overset{n}{\sum}}}{\displaystyle {\int}_{0}^{\infty}}\left\{{X}_{i}\left(t\right)-\frac{{\displaystyle {\sum}_{i=1}^{n}}\text{\hspace{0.05em}}{{X}^{\prime}}_{j}\left(t\right){Y}_{j}\left(t\right)\mathrm{exp}\left({\vartheta}^{\prime}\left(t\right){X}_{j}\left(t\right)\right)}{{\displaystyle {\sum}_{i=1}^{n}}\text{\hspace{0.05em}}{Y}_{j}\left(t\right)\mathrm{exp}\left({\vartheta}^{\prime}\left(t\right){X}_{j}\left(t\right)\right)}\right\}\text{d}{N}_{i}\left(t\right).$ (15)

Someone may use computationally easiest alternative versions of (12) which were first mentioned by [5] and lately by [11] .

Finally, the survival function of T given possibly time-varying covariates $X\left(t\right)$ can easily be derived from the model (5) as follows.

${S}_{T\mathrm{|}X\left(t\right)}\left\{t\mathrm{|}x\left(t\right)\right\}=\mathrm{exp}\left\{-\Phi \left({\displaystyle {\int}_{0}^{t}}\text{\hspace{0.05em}}\text{d}{\Lambda}_{\epsilon}\left\{\mathcal{H}\left(u\right)+{\vartheta}^{\prime}\left(u\right)X\left(u\right)\right\}\right)\right\}\mathrm{.}$ (16)

Therefore, the cumulative hazard function is given by

${\stackrel{^}{\Lambda}}_{T}\left(t\mathrm{,}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\mathrm{,}\stackrel{^}{\mathcal{H}}\left(\mathrm{.;}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\right)\mathrm{|}\left({X}_{i}\left(t\right)\right)\right)=-\mathrm{log}{\stackrel{^}{S}}_{T}\left(t\mathrm{,}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\mathrm{,}\stackrel{^}{\mathcal{H}}\left(\mathrm{.;}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\right)\mathrm{|}\left({X}_{i}\left(t\right)\right)\right)\mathrm{,}$ (17)

thus,

$\begin{array}{l}{\stackrel{^}{\Lambda}}_{T}\left(t\mathrm{,}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\mathrm{,}\stackrel{^}{\mathcal{H}}\left(\mathrm{.;}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\right)\mathrm{|}\left({X}_{i}\left(t\right)\right)\right)\\ =-\mathrm{log}\left[\mathrm{exp}\left\{-\Phi \left({\displaystyle {\int}_{0}^{t}}\text{\hspace{0.05em}}\text{d}{\Lambda}_{\epsilon}\left\{\stackrel{^}{\mathcal{H}}\left(u\mathrm{;}{\stackrel{^}{\vartheta}}_{n}\left(u\right)\right)+{{\stackrel{^}{\vartheta}}^{\prime}}_{n}\left(u\right){X}_{i}\left(u\right)\right\}\right)\right\}\right]\mathrm{.}\end{array}$ (18)

Thus, the true induced intensity (hazard) function for failure time T given possibly time-varying covariates $X\left(t\right)$ is the derivative of the true cumulative intensity function of Equation (18) which is defined as

${\lambda}_{T}\left(t\mathrm{,}\vartheta \left(t\right)\mathrm{,}\mathcal{H}\left(\mathrm{.;}\vartheta \left(t\right)\right)\mathrm{|}\left({X}_{i}\left(t\right)\right)\right)=\frac{\text{d}}{\text{d}t}{\Lambda}_{T}\left(t\mathrm{,}\vartheta \left(t\right)\mathrm{,}\mathcal{H}\left(\mathrm{.;}\vartheta \left(t\right)\right)\mathrm{|}\left({X}_{i}\left(t\right)\right)\right),$ (19)

therefore, to ease the notations without lose of truth, here we propose some representations

${\lambda}_{T}\left(t\mathrm{,}\vartheta \left(t\right)\mathrm{,}\mathcal{H}\left(\mathrm{.;}\vartheta \left(t\right)\right)\mathrm{|}\left({X}_{i}\left(t\right)\right)\right)=\stackrel{\u02dc}{h}\left(t\mathrm{,}\vartheta \mathrm{,}\mathcal{H}\mathrm{|}\left({X}_{i}\left(t\right)\right)\right)\ast \stackrel{\dot{}}{\mathcal{H}}\left(t\right)\mathrm{,}\stackrel{\dot{}}{\mathcal{H}}\left(t\right)=\frac{\partial}{\partial t}\mathcal{H}\left(t\right)\mathrm{,}$ (20)

where

$\begin{array}{l}\stackrel{\u02dc}{h}\left(t\mathrm{,}\vartheta \left(t\right)\mathrm{,}\mathcal{H}\mathrm{|}{X}_{i}\left(t\right)\right)\\ =\stackrel{\dot{}}{\Phi}\left({\displaystyle {\int}_{0}^{t}}\text{\hspace{0.05em}}\text{d}{\Lambda}_{\epsilon}\left\{\stackrel{^}{\mathcal{H}}\left(u\mathrm{;}{\stackrel{^}{\vartheta}}_{n}\left(u\right)\right)+{{\stackrel{^}{\vartheta}}^{\prime}}_{n}\left(u\right){X}_{i}\left(u\right)\right\}\right)\ast {\lambda}_{\epsilon}\left\{\stackrel{^}{\mathcal{H}}\left(u\mathrm{;}{\stackrel{^}{\vartheta}}_{n}\left(u\right)\right)+{{\stackrel{^}{\vartheta}}^{\prime}}_{n}\left(u\right){X}_{i}\left(u\right)\right\}\mathrm{,}\end{array}$ (21)

where $\stackrel{\dot{}}{\Phi}\left(t\right)=\frac{\partial}{\partial t}\Phi \left(t\right)$ .

Now, we set a zero-mean martingale process with respective filtration $\mathcal{F}$ of complete $\sigma \left\{{Y}_{i}\left(u\right)\mathrm{,}{N}_{i}\left(u\right)\mathrm{,}X\left(u\right)\mathrm{,0}\le u\le t\right\}$ as

${M}_{i}\left(t\right)={N}_{i}\left(t\right)-{\displaystyle {\int}_{0}^{t}}\text{\hspace{0.05em}}{Y}_{i}\left(u\right){\lambda}_{T}\left(t\mathrm{,}\vartheta \left(t\right)\mathrm{,}\mathcal{H}\left(\mathrm{.;}\vartheta \left(t\right)\right)\mathrm{|}\left({X}_{i}\left(t\right)\right)\right)\text{d}u\mathrm{.}$ (22)

Thus, by imposing at Lemma 1, we modify the estimating Equation (12) and Equation (13) as

$\begin{array}{l}{Q}_{1}\left(\vartheta \left(t\right)\mathrm{,}\mathcal{H}\left(\mathrm{.}\right)\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum}}}{\displaystyle {\int}_{0}^{t}}{X}_{i}\left(u\right)\left[\text{d}{N}_{i}\left(u\right)-{Y}_{i}\left(u\right){\lambda}_{T}\left(u\mathrm{,}\vartheta \left(u\right)\mathrm{,}\mathcal{H}\left(\mathrm{.;}\vartheta \left(u\right)\right)\mathrm{|}\left({X}_{i}\left(u\right)\right)\right)\text{d}u\right]=\mathrm{0,}\end{array}$ (23)

$\begin{array}{l}{Q}_{2}\left(\vartheta \left(t\right)\mathrm{,}\mathcal{H}\left(\mathrm{.}\right)\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum}}}\left[\text{d}{N}_{i}\left(t\right)-{Y}_{i}\left(t\right){\lambda}_{T}\left(t\mathrm{,}\vartheta \left(t\right)\mathrm{,}\mathcal{H}\left(\mathrm{.;}\vartheta \left(t\right)\right)\mathrm{|}{X}_{i}\left(t\right)\right)\text{d}t\right]=0\left(t\ge 0\right)\mathrm{.}\end{array}$ (24)

3. Large Sample Theory and Conditions

Some regularity conditions are necessarily imposed here.

C1: The covariate vectors are bounded in the sense that $Pr\left(\Vert \text{\hspace{0.05em}}\mathrm{.}\text{\hspace{0.05em}}\Vert <L\right)=1$ for some constant $L>0$ and the possibly time-varying covariate $X\left(t\right)$ has a uniformly bounded variation on $\left[\mathrm{0,}\mathcal{T}\right]$ and its left limit exists with any t where $\mathcal{T}$ is the maximum follow-up time.

C2: The true value of $\vartheta \left(t\right)$ denote by ${\vartheta}_{0}\left(t\right)$ , lies in the interior of a known compact set $\Theta $ in ${R}^{k}$ and the true value of $\mathcal{H}\left(\mathrm{.}\right)$ denote by ${\mathcal{H}}_{0}\left(\mathrm{.}\right)$ is continuously positive differentiable on the closed interval $\left[\mathrm{0,}\mathcal{T}\right]$ .

C3: The transformation
$\Phi \left(\mathrm{.}\right)$ is at least thrice continuously differentiable on interval
$\left[\mathrm{0,}\infty \right)$ with
$\Phi \left(0\right)=0$ and
${\Phi}^{k}(.)>0,\Phi \left(\infty \right)=\infty $ , and
${\mathrm{sup}}_{u\ge 0}\left|{\Phi}^{k}(.)\right|<\infty ,\mathrm{}k=1,2,3$ , where
${\Phi}^{k}(.)$ denotes k^{th} derivatives of
$\Phi \left(\mathrm{.}\right)$ .

C4: ${\Lambda}_{\epsilon}\left(\mathrm{.}\right)$ is a strictly increasing positive function on interval $\left[\mathrm{0,}\mathcal{T}\right]$ and ${\Lambda}_{\epsilon}\left(\mathrm{.}\right)$ is continuously differentiable.

C5: For any given finite scalar k, ${\lambda}_{\epsilon}\left(\mathrm{.}\right)$ is strictly positive and ${\stackrel{\dot{}}{\lambda}}_{\epsilon}\left(\mathrm{.}\right)$ is bounded and continuously differentiable on interval $\left(-\infty \mathrm{,}k\right]$ , where the superscript dot always refers derivatives.

C6: Both the variance covariance matrices $\Psi $ and $\stackrel{\u02dc}{\Sigma}$ are nonsingular.

Theorem 1: Under some suitable regularity conditions C1-C6 in order to ensure CLT for counting process martingale holds, ${\stackrel{^}{\vartheta}}_{n}\left(t\right)$ is consistent estimator of the true parameter ${\vartheta}_{0}\left(t\right)$ , i.e.;

${n}^{\frac{1}{2}}\left({\stackrel{^}{\vartheta}}_{n}\left(t\right)-{\vartheta}_{0}\left(t\right)\right)\stackrel{d}{\to}N\left(\mathrm{0,}\left({\Psi}^{-1}\right)\stackrel{\u02dc}{\Sigma}{\left({\Psi}^{-1}\right)}^{\prime}\right)\mathrm{.}$ (25)

Thus, similar to [5] [9] [12] and others, the asymptotic variance of estimator ${\stackrel{^}{\vartheta}}_{n}\left(t\right)$ can be estimated consistently by estimating $\Psi $ and $\stackrel{\u02dc}{\Sigma}$ consistently.

Theorem 2: Under some suitable conditions C1-C6, ${\stackrel{^}{\mathcal{H}}}_{n}\left(t\mathrm{;}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\right)$ is consistent under the metric $d\left(\cdot \mathrm{,}\cdot \right)$ , where for any two nondecreasing functions ${\stackrel{^}{\mathcal{H}}}_{n}\left(t\mathrm{;}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\right)$ and $\mathcal{H}\left(t\right)$ on interval $\left[\mathrm{0,}\infty \right)$ such that ${\stackrel{^}{\mathcal{H}}}_{n}\left(0\right)={\mathcal{H}}_{0}\left(0\right)=-\infty $ ,

$d\left({\stackrel{^}{\mathcal{H}}}_{n}\left(t\mathrm{;}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\right)\mathrm{,}{\mathcal{H}}_{0}\left(t\right)\right)=\mathrm{sup}\left(\left|{\stackrel{^}{\mathcal{H}}}_{n}\left(t\mathrm{;}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\right)-{\mathcal{H}}_{0}\left(t\right)\right|\mathrm{:}t\in \left[b\mathrm{,}\mathcal{T}\right]\right)\stackrel{P}{\to}\mathrm{0,}$ (26)

for any fixed $b\in \left(\mathrm{0,}\mathcal{T}\right]$ , where $\mathcal{T}=\mathrm{inf}\left(t\mathrm{:}P\left(\stackrel{\u02dc}{Y}>t\right)=0\right)$ . The following theorem establishes the asymptotic distribution of the estimated distribution.

4. Simulation Study

The data is simulated from Cox model in four cases, such as with fixed covariates, with fixed covariates and time-varying coefficients, with time-varying covariates and with time-varying covariates and time-dependent effects simultaneously. The data was generated by using sim.survdata() under R package “coxed” based on the flexible hazard methods described by [13] . The survival time data with three continous covariates was generated with sample size n = 200 and maximum duration 50 units using sim.survdata(). By default sim.survdata() generates the survival time and three covariates from standard normal distribution. However, we can adjust for other characteristics of covariates from different distributions for fixed covariate case.

Required data structure for time-dependent covariate is technically different from the survival data structure with baseline covariates. The dependent variable for Cox model in survival data can be arranged by using “Surv()” function in survival package of R software. Commonly it has two arguments survival time and a censoring time variables. However, for in the case of time-varying covariates the survival time variable setup is divided in to two sections referring start and end of discrete intervals, which in turn permits a covariate to be measured in different values across different intervals for the same observations. Thus, in the case of time dependent covariates, we set type = “tvc” in “sim.survdata()” function to generated survival time data with time varying covariates. Then the survival durations are generated again using proportional hazards, and are passed to the “permalgorithm()” functionin the “permAlgo” package to generate the time-varying data structure [14] . In the case of time-dependent covariates, the type = “tvc” option of sim.survdata does not allow to use user supplied data for the covariates, as a time-varying covariate is expressed overtime frames which themselves convey part of the variation of the times, and then the time is generated [15] .

The usual proportionality assumptions of Cox proportional hazard model fails when the coefficient effect varies through over time. The data for time-dependent coefficients can similarly generated using sim.survdata() function by setting the type = “tvbeta” option inside the function. Whenever this option sets, the first coefficient, whether coefficients are user-supplied or randomly generated, is interacted with natural logarithm of the time counter from 1 to maximum time T [15] . Then the sim.survdata() function generates survival time from proportional hazards model, and saves the coefficients in designed matrix form to allow their dependence on time. So to generate the data with the time-dependent coefficients we set type = “tvbeta”.

The data for more flexible and general cox model with the time-dependent coefficients and the time-dependent covariates can similarly generated using sim.survdata() function by setting the type = c(“tv”, ”tvbeta”) option inside the function. Finally, semiparametric transformation models are applied for the simulated data. The different models were compared based on their performance in precision.

4.1. Computational Algorithm

Since we have more than one unknown items to be estimated, it is necessary to apply some sophisticated iterative algorithms to handle the iteration problem. Thus, in this paper expectation-maximization (EM) algorithm is proposed to estimate unknown true parameter ${\vartheta}_{0}\left(t\right)$ and nondecreasing monotone function ${\mathcal{H}}_{0}\left(\mathrm{.;}\vartheta \left(t\right)\right)$ . In this concept, it is necessary to fix one of them and estimate the another one and in terms of the fixed one and vice versa. Therefore, as it was done in [5] , it is not difficult to show the unique solution of (12), (13) in H, for every fixed value of $\vartheta \left(t\right)$ . Consequently, Equation (3) and Equation (5) logically suggest the following iterative algorithms for computing $\left({\stackrel{^}{\vartheta}}_{n}\left(t\right)\mathrm{,}{\stackrel{^}{\mathcal{H}}}_{n}\left(\mathrm{.;}{\stackrel{^}{\vartheta}}_{n}\left(t\right)\right)\right)$ .

Step 0: Opt an initial value of $\vartheta $ , denoted by ${\vartheta}^{\left(0\right)}$ .

Step 1: For each ${t}_{k}$ , obtain ${\vartheta}^{\left(0\right)}\left(t\right)$ and ${\mathcal{H}}^{\left(0\right)}\left({t}_{1}\right)$ by solving Equation (12) and Equation (13) by setting

$\underset{i=1}{\overset{n}{\sum}}}\text{\hspace{0.05em}}{Y}_{i}\left({t}_{1}\right){\lambda}_{\epsilon}\left\{{\vartheta}^{\prime}\left(t\right){X}_{i}\left({t}_{1}\right)+\mathcal{H}\left({t}_{1}\right)\right\}\stackrel{\dot{}}{\Phi}\left({\displaystyle {\int}_{0}^{t}}\left[{\Lambda}_{\epsilon}\left\{{\vartheta}^{\prime}\left(u\right){X}_{i}\left({u}_{1}\right)+\mathcal{H}\left({u}_{1}\right)\right\}\right]\right)=\mathrm{1,$

with $\vartheta \left(t\right)={\vartheta}^{0}\left(t\right)$ . Then obtain ${\mathcal{H}}^{\left(0\right)}\left({t}_{1}\right)$ , for $k=2,3,\cdots ,K$ , one-by-one by solving the equation

$\begin{array}{l}{\displaystyle \underset{i=1}{\overset{n}{\sum}}}\text{\hspace{0.05em}}{Y}_{i}\left({t}_{1}\right){\lambda}_{\epsilon}\left\{{\vartheta}^{\prime}\left(t\right){X}_{i}\left({t}_{1}\right)+\mathcal{H}\left({t}_{1}\right)\right\}\stackrel{\dot{}}{\Phi}\left({\displaystyle {\int}_{0}^{t}}\left[{\Lambda}_{\epsilon}\left\{{\vartheta}^{\prime}\left(u\right){X}_{i}\left({u}_{1}\right)+\mathcal{H}\left({u}_{1}\right)\right\}\right]\right)\\ =1+{\displaystyle \underset{i=1}{\overset{n}{\sum}}}\text{\hspace{0.05em}}{Y}_{i}\left({t}_{k-1}\right){\lambda}_{\epsilon}\left\{{\vartheta}^{\prime}\left(t\right){X}_{i}\left({t}_{k-1}\right)+\mathcal{H}\left({t}_{k-1}\right)\right\}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\cdot \stackrel{\dot{}}{\Phi}\left({\displaystyle {\int}_{0}^{t}}\left[{\Lambda}_{\epsilon}\left\{{\vartheta}^{\prime}\left(t\right){X}_{i}\left({u}_{k-1}\right)+\mathcal{H}\left({u}_{k-1}\right)\right\}\right]\right).\end{array}$

Step 2: Then obtain new estimate of $\vartheta \left(t\right)$ by solving (12) with $\stackrel{^}{\mathcal{H}}\left({t}_{k}\right)={\mathcal{H}}^{\left(0\right)}\left({t}_{k}\right)$ as obtained in Step 1.

Step 3: Set ${\vartheta}^{\left(0\right)}\left(t\right)$ to be the estimator obtained in Step 2 and repeat Steps 1 and 2 until prescribed convergence criteria are met based on Equation (12) and Equation (13).

4.2. Numerical Results

This subsection explores the numerical results based on simulation studies through figures and numerical analysis. This numerical result is expected to evaluate the performance of the proposed model.

Figure 1 illustrates about the baseline characteristics of survival data. The top panel of the figure refers the feature of probability density function, cumulative distribution function, hazard function, and cumulative hazard function of failure time. The bottom panel shows the feature of simulated duration in terms of histogram of failure time or duration, linear predictor and exponentiated linear predictors respectively. The left panel of Figure 1 is when the survival data are

(a) (b)

Figure 1. Plots of baseline feature of simulated survival data. (a) Plot with 25% censoring rate; (b) Plot with 45% censoring rate.

assumed to have 25% censoring rate and the right side panel is when the survival data are assumed to have 45% censoring rate.

Table 1 illustrates the results of simulation based on four different cases under special cases of semiparametric transformation models. The result has shown, the performance of the model reduces as censoring rate increases. The standard errors in the bracket indicated the precision level of the estimators. The estimators with small standard errors have high precision. In these simulations, the effect of time-varying coefficient did not improve the model performance. However, the effect of time-varying covariates did improve the performance of the model.

5. Conclusions

The study is basically concerned on comparisons of the semiparametric transformation models with and without the effect of the time on covariates and coefficients. The summary review of other works was done and the result of simulation was included to come up with reasonable review of the study. The data were generated in four different cases under the “sim.survdata()” function of R package called “coxed”. Then the results of semiparametric transformation models for four types of simulation studied were compared based. Three special cases of semiparametric transformation models such as PH, PO and model when r = 0.5 were considered.

The results have shown that the semiparametric transformation models with time-dependent covariates did relatively better perform with small standard errors. However, the effect of time-varying coefficient did not improve the performance of the semiparametric transformation models in our simulation studies. The last

Table 1. Estimates of Regression Coefficients with their respective standard errors in the brackets for Semiparametric Transformation models for n = 200. TCV and TVbeta refers time-varying covariates and time-varying coefficients.

two cases such as the semiparametric transformation models with time-varying covariates and both time-varying covariates and time-varying coefficients have shown better performance. Therefore, we can give the general conclusion that when the proportionality assumption fails to fulfill, incorporating the time-varying coefficient effect in the model is advisable. Considering only baseline covariate may not be always true; because there is the time when the covariate changes throughout the time. Therefore, incorporating time-varying covariate in the model may help us to get reasonable results. Sometimes it can be happened that both covariate and coefficient effect changes over time. Thus, incorporating both time-varying covariates and time-varying coefficients shall give us more reasonable results.

References

[1] Wang, X. and Zhou, X.H. (2017) Semiparametric Maximum Likelihood Estimation for the Cox model with Length-Biased Survival Data. Journal of Statistical Planning and Inference, 196, 163-173.

https://doi.org/10.1016/j.jspi.2017.11.004

[2] Cox, D.R. (1972) Regression Models and Life-Tables. In: Kotz, S. and Johnson, N.L., Eds., Springer, New York, NY.

[3] Horowitz, J.L. (1998) Semiparametric Methods in Economics. In: Kotz, S. and Johnson, N.L., Eds., Lecture Notes in Statistics, Springer, New York.

[4] Cheng, S.C., Wei, L.J. and Ying, Z. (1995) Analysis of Transformation Models with Censored Data. Biometrika, 82, 835-845.

https://doi.org/10.1093/biomet/82.4.835

[5] Chen, K.N., Jin, Z.Z. and Ying, Z.L. (2002) Semiparametric Analysis of Transformation Models with Censored Data. Biometrika, 89, 659-668.

https://doi.org/10.1093/biomet/89.3.659

[6] Zeng, D.L. and Lin, D.Y. (2006) Efficient Estimation of Semiparametric Transformation Models for Counting Processes. Biometrika, 93, 627-640.

https://doi.org/10.1093/biomet/93.3.627

[7] Chen, L., Lin, D.Y. and Zeng, D.L. (2012) Checking Semiparametric Transformation Models with Censored Data. Biostatistics, 13, 18-31.

https://doi.org/10.1093/biostatistics/kxr017

[8] Zeng, D. and Lin, D.Y. (2007) Maximum Likelihood Estimation in Semiparametric Regression Models with Censored Data. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69, 507-564.

https://doi.org/10.1111/j.1369-7412.2007.00606.x

[9] Kim, J.P., Lu, W.B., Sit, T. and Ying, Z.L. (2013) A Unified Approach to Semiparametric Transformation Models Under General Biased Sampling Schemes. Journal of the American Statistical Association, 108, 217-227.

https://doi.org/10.1080/01621459.2012.746073

[10] Sinha, S. and Ma, Y.Y. (2014) Semiparametric Analysis of Linear Transformation Models with Covariate Measurement Errors. Biometrics, 70, 21-32.

https://doi.org/10.1111/biom.12119

[11] Ma, H.J., Qiu, Z.P. and Zhou, Y. (2016) Semiparametric Transformation Models with Length-Biased and Right-Censored Data under the Case-Cohort Design. Statistics and Its Interface, 9, 213-222.

https://doi.org/10.4310/SII.2016.v9.n2.a8

[12] Shen, P.S. (2017) Semiparametric Analysis of Transformation Models with Dependently Lefttruncated and Right-Censored Data. Communications in Statistics— Simulation and Computation, 46, 2474-2487.

https://doi.org/10.1080/03610918.2015.1048879

[13] Harden, J.J. and Kropko, J. (2018) Simulating Duration Data for the Cox Model. Political Science research and Methods, 1-8.

https://doi.org/10.1017/psrm.2018.19

[14] Kropko, J. and Harden, J.J. (2018) How to Simulate Survival Data with the sim.survdata Function.

https://cran.r-project.org/web/packages/coxed/vignettes/simulating_survival_data.html

[15] Sylvestre, M.P. and Abrahamowicz, M. (2018) Comparison of Algorithms to Generate Event Times Conditional on Time-Dependent Covariates. Statistics in Medicine, 27, 2618-2634.

https://doi.org/10.1002/sim.3092