Regression Analysis of Multivariate Interval-Censored Failure Time Data under Transformation Model with Informative Censoring

Abstract

We consider a regression analysis of multivariate interval-censored failure time data where the censoring may be informative. To address this, an approximated maximum likelihood estimation approach is proposed under a general class of semiparametric transformation models, and in the method, the frailty approach is employed to characterize the informative interval censoring. For the implementation of the proposed method, we develop a novel EM algorithm and show that the resulting estimators of the regression parameters are consistent and asymptotically normal. To evaluate the empirical performance of the proposed estimation procedure, we conduct a simulation study, and the results indicate that it performs well for the situations considered. In addition, we apply the proposed approach to a set of real data arising from an AIDS study.

1. Introduction

In this paper, we consider a regression analysis of multivariate interval-censored failure time data where the censoring may be informative. Interval-censored data arise when the failure time of interest is known or observed only to belong to an interval instead of being observed exactly (Finkelstein, 1986 [1]; Sun, 2006 [2]). It is apparent that one can treat right-censored data as a special case of interval-censored data. Multivariate interval-censored data occur if a failure time study involves more than one related failure time of interest for which only interval-censored data are available. Among others, one can often face multivariate interval-censored data in epidemiological studies and clinical trials. Informative censoring occurs if the censoring mechanism or the underlying process generating observations is related to the failure times of interest (Kalbfleisch and Prentice, 2002 [3]; Sun, 2006 [2]).

An example of informative censoring is given by a clinical trial or periodic study on a failure event such as death for which some symptoms may occur before the occurrence of the event. For the situation, the study subject may tend to pay more clinical visits when the symptoms occur rather than following the pre-specified schedule. Many authors have pointed out’ that with informative censoring, the analysis that ignores it could lead to serious biased estimators or analysis results (Wang et al., 2010 [4]; Sun, 2006 [2]; Zhang et al., 2005 [5], 2007 [6]). For example, Sun (1999) [7] studied the issue for univariate current status data, a special case of interval-censored data where the observed interval includes either zero or infinity, and showed that the analysis could yield misleading results if the informative censoring is ignored or treated to be non-informative censoring.

A large amount of literature has been established for the regression analysis of multivariate interval-censored failure time data or their special cases, multivariate current status data and bivariate interval-censored data (Chen et al., 2007 [8], 2009 [9], 2013 [10]; Goggins and Finkelstein, 2000 [11]; Shen, 2015 [12]; Tong et al., 2008 [13]; Wang et al., 2008 [14]; Zeng et al., 2017 [15]; Zhang et al., 2009 [16]; Zhou et al., 2017 [17]). For this, three types of methods are commonly used, including the copula model approach, the marginal model-based approach and the frailty model-based method. The first employs various copula models to characterize the relationship between correlated failure times of interest (Wang et al., 2008 [14]; Zhang et al., 2009 [16]), and among others, Sun and Ding (2019) [18] discussed this for bivariate cases under the framework of the two-parameter Archimedean copula model. The second mainly focuses on the marginal distribution and puts no assumption on the correlation between the failure times of interest (Wei et al. 1989 [19]). The authors who developed such methods include Chen et al. (2007) [8], Chen et al. (2013) [10] and Tong et al. (2008) [13].

The frailty model-based approach generally employs the frailty or latent variable to model the correlation between the correlated failure times. It has the advantage of allowing one to directly estimate the correlation. One main shortcoming of most of the existing methods for multivariate interval-censored data is that they assume independent or non-informative interval censoring, and it is apparent that this may not hold in practice as discussed above. In this paper, we will adopt the frailty model-based approach to develop a new estimation procedure that allows for dependent or informative interval censoring.

Several authors have considered a regression analysis of univariate informatively interval-censored failure time data. For example, Zhang (2005) [5], Wang et al. (2010) [4] and Wang et al. (2018) [20] investigated the problem for current status data, case II interval-censored data and case K interval-censored data, respectively. Case II means that each study subject is observed only twice, while case K refers to the situation where each subject is observed at a sequence of observation times, which is much more general than others (Sun, 2006) [2]. As mentioned above, most of the existing methods for multivariate interval-censored data apply only to the situation with independent interval censoring, except Yu et al. (2022) [21]. Yu et al. (2022) [21] only considered case II interval-censored data under the additive hazards model. In this paper, the focus will be on case K multivariate interval-censored data with informative censoring and the proposed methods apply to much more general situations than Yu et al. (2022) [21].

More specifically, in Section 2, some notation and assumptions will be first introduced as well as the data structures. In the proposed method, we will focus on the case where the failure time of interest marginally follows a general class of semiparametric transformation models. The proposed approximated maximum likelihood estimation procedure will be presented in Section 3, and for the implementation of the proposed method, a novel EM algorithm will be developed. The asymptotic properties of the resulting estimators of the regression parameters will be given in Section 4. Section 5 will present some simulation results obtained from a study performed to evaluate the performance of the proposed method, and they indicate that it performs as expected. In Section 6, we apply the proposed methodology to a set of real data arising from an AIDS clinical trial, and Section 7 contains some discussion and concluding remarks.

2. Assumptions and Background

In this section, we first introduce some notation and background and then describe the model and data structure. Suppose that there is a failure time study consisting of n independent subjects and concerning M failure events of interest that may be related. Define $T_{im}$ to be the failure time of interest and $X_{im}$ a p-dimensional vector of covariates both related to the ith subject and the event m. Furthermore, for each subject, suppose that there exists a sequence of potential observation times $U_{i0}=0<U_{i1}<\dots <U_{i{K}_i^{\ast }}$ and a follow-up or stopping time ${\tau }_i$, where $K_i^{\ast }$ denotes the number of potential observations, $i=1,\dots ,n$. For simplicity, we assume that for each subject, the observation times for different failure events are the same and the proposed method below can be easily generalized to more general situations.

For subject i, define the point process ${\tilde{N}}_i\left(t\right)={\sum }_{j=1}^{K_i^{\ast }}I(U_{ij}\le t)$, describing the observation process on the subject that jumps only at the observation times, $i=1,\dots ,n$. Note that for the situation considered here, we have $M+1$ processes, the M underlying failure time processes of interest and the observation process ${\tilde{N}}_i\left(t\right)$, and as mentioned above, the focus below will be on the case where they may be related (Ma et al., 2015 [22]; Wang et al., 2016 [23]; Zhang et al., 2007 [6]). To describe their relationships and the possible covariate effects on them, we assume that there exists a vector of latent variables $b_i$ and another latent variable $u_i$ with mean zero, and given $X_{im}$, $Z_{im}$, $b_i$ and $u_i$, the cumulative hazard function of $T_{im}$ has the form

$${\Lambda }_{im}\left(t\right|X_{im},Z_{im},b_i,u_i)=G_m\{exp(X_{im}^T{\beta }_{xm}+Z_{im}^T{b}_i+u_i{\beta }_{um}){\Lambda }_m\left(t\right)\}.$$

(1)

Here, $G_m(.)$ is a known strictly increasing transformation function, ${\Lambda }_m(.)$ is an unknown baseline cumulative hazard function, $Z_{im}$ contains 1 and part of the covariates $X_{im}$, and ${\beta }_m={({\beta }_{xm}^T,{\beta }_{um})}^T$ denotes the vector of unknown regression parameters.

For the observation process, it will be assumed that ${\tilde{N}}_i\left(t\right)$ is a non-homogeneous Poisson process satisfying

$${\lambda }_{ih}\left(t\right|X_i,u_i)={\lambda }_{0h}\left(t\right)exp\left(\sum _{m=1}^M{X}_{im}^T{\alpha }_m+u_i\right)={\lambda }_{0h}\left(t\right)exp\left(X_i^T\alpha +u_i\right)$$

(2)

for the intensity function given $X_{im}$ and $u_i$. Here, ${\lambda }_{0h}$ is an unknown continuous baseline intensity function, $X_i={(X_{i1}^T,\dots ,X_{iM}^T)}^T$, and ${\alpha }^T=({\alpha }_1^T,\dots ,{\alpha }_M^T)$, which is a vector of regression parameters as ${\beta }_m$. In the following, it will be assumed that given $X_{im}$, $b_i$ and $u_i$, $T_{i1},\dots ,T_{iM}$ are independent, and given $X_{im}$ and $u_i$, $T_i$ and ${\tilde{N}}_i$ are independent. Moreover, ${\tau }_i$ is independent of $T_i$ and ${\tilde{N}}_i$. We point out that models (1) and (2) with $u_i=0$ have been commonly used in the analysis of failure-time data (Klein and Moescherger, 2003 [24]) and the analysis of event history data (Cook and Lawless, 2007 [25]), respectively. The parameter ${\beta }_{um}$ denotes the degree of the correlation between the failure-time process and the observation process, and they will be independent if ${\beta }_{um}=0$.

The semiparametric transformation model (1) with $b_i=0$ and $u_i=0$ is quite general and can give many specific models. In particular, one can express it as a class of frailty-induced transformations

$$G_m\left(x\right)=-log{\int }_0^{\infty }exp(-xt)f_m\left(t\right)dt.$$

In the above, $f_m\left(t\right)$ denotes the density function of a frailty variable with support $[0,\infty ]$. By setting $f_m\left(t\right)$ to be the gamma density with mean 1 and variance $r_m$, it gives the class of logarithmic transformations $G_m\left(x\right)=r_m^{-1}log(1+r_{m}x)$ with $r_m>0$ (Chen et al., 2002 [26]). In particular, it yields the proportional odds model with $r_m=1$ or $G_m\left(x\right)=log(1+x)$ and gives the proportional hazards model with $r_k=0$ or $G_m\left(x\right)=x$.

To describe the observed data, define ${\delta }_{imj}=I(U_{ij-1}<T_{im}\le U_{ij})$, indicating if the failure time $T_{im}$ belongs to the interval $(U_{ij-1},U_{ij}]$. In the following, it will be assumed that the observed data have the form

$$\mathcal{O}=\{O_i=({\tau }_i,X_{im},Z_{im},U_{ij},{\delta }_{imj},m=1,\dots ,M,j=1,\dots ,K_i^{\ast });i=1,\dots ,n\},$$

where $K_i^{\ast }={\tilde{N}}_i\left({\tau }_i\right)$. That is, we observe case K interval-censored data (Sun, 2006).

3. Maximum Likelihood Estimation

3.1. Estimation Procedure

Now, we discuss inference about models (1) and (2), and for this, we will propose a two-step or an approximate maximum likelihood estimation procedure by following Huang and Wang (2004) [27] and Wang et al. (2016) [23]. More specifically, we will first consider the estimation of model (2) and then estimation of ${\varphi }^T=({\beta }_1^T,\dots ,{\beta }_M^T)$, the parameter of interest. The first step will be based on the following two facts. One is that one can easily show that $K_i^{\ast }$ follows the Poisson distribution with the mean

$${\Lambda }_{ih}({\tau }_i;X_i,u_i)={\Lambda }_{0h}\left({\tau }_i\right)exp(X_i^T\alpha +u_i)$$

given $X_i$ and $u_i$, $i=1,\dots ,n$. The other is that the observation times $U_{i1},\dots ,U_{i{K}_i^{\ast }}$ can be seen as the order statistics of a set of i.i.d. random variables with the density function

$$\begin{array}{c}\hfill \pi \left(t\right)=\frac{{\lambda }_{0h}\left(t\right)exp(X_i^T\alpha +u_i)}{{\Lambda }_{0h}\left({\tau }_i\right)exp(X_i^T\alpha +u_i)}I(0\le t\le {\tau }_i)=\frac{{\lambda }_{0h}\left(t\right)}{{\Lambda }_{0h}\left({\tau }_i\right)}I(0\le t\le {\tau }_i).\end{array}$$

One can see that the function above does not depend on neither $X_i$ nor $u_i$, which suggests that the function ${\Lambda }_{0h}\left(t\right)$ can be estimated by

$${\widehat{\Lambda }}_{0h}\left(t\right)=\prod _{s\left(l\right)>t}(1-\frac{d_{\left(l\right)}}{R_{\left(l\right)}}).$$

In the above, the $s_{\left(l\right)}$’s denote the ordered and distinct values of observation times $\left\{U_{ik}\right\}$, $d_{\left(l\right)}$ the number of the observation times equal to $s_{\left(l\right)}$, and $R_{\left(l\right)}$ the number of observation times satisfying $U_{ik}\le s_{\left(l\right)}\le {\tau }_i$ among all subjects.

Under the assumptions above, it is easy to show that $E[K_i^{\ast };X_{i.},u_i,{\tau }_i]={\Lambda }_{0h}\left({\tau }_i\right)exp$$(X^T\alpha +u_i)$. This yields

$$E_{u_i}\left[E[K_i^{\ast }{\Lambda }_{0h}^{-1}\left({\tau }_i\right);X_{im},u_i,{\tau }_i]\right]=E\left(e^{u_i}\right)exp\left(X_i^T\alpha \right),$$

and a class of estimating equations

$$\sum _{i=1}^n{\omega }_i{X}_i(K_i^{\ast }{\widehat{\Lambda }}_{0h}^{-1}\left({\tau }_i\right)-E\left(e^{u_i}\right)exp\left(X_i^T\alpha \right))=0$$

for estimation of ${\alpha }_m$, $m=1,\dots ,M$ with the ${\omega }_i$’s being some weights. Let ${\widehat{\alpha }}_m$ denote the estimator of ${\alpha }_m$ given by the estimating equations above, which suggests that one can naturally estimate $u_i$ by

$${\widehat{u}}_i=log\left\{\frac{K_i^{\ast }}{{\widehat{\Lambda }}_{0h}\left({\tau }_i\right)exp\left(X_i^T\alpha \right)}\right\}.$$

Now, consider estimation of $\varphi$ as well as model (1). For this, note that if the $u_i$’s were known, it would be natural to maximize the likelihood function

$$\begin{array}{ccc}\hfill L_n(\varphi ,\Lambda ,\gamma |u_i^{\prime }s)=& & \prod _{i=1}^n\int \prod _{m=1}^M\prod _{j=1}^{K_i^{\ast }}\left\{exp(-G_m[{\int }_0^{U_{i,j-1}}exp\{x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i\}d{\Lambda }_m\left(t\right)])-\right.\hfill \\ & & {\left.exp(-G_m[{\int }_0^{U_{ij}}exp\{x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i\}d{\Lambda }_m\left(t\right)])\right\}}^{{\Delta }_{imj}}\hfill \\ & & exp{(-G_m[{\int }_0^{U_{i{K}_i^{\ast }}}exp\{x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i\}d{\Lambda }_m\left(t\right)])}^{1-{\sum }_{j=1}^{K_i^{\ast }}{\Delta }_{imj}}{f}_b\left(b_i\right|\gamma )d{b}_i,\hfill \end{array}$$

where $\Lambda =({\Lambda }_1,\dots ,{\Lambda }_M)$, $x_{im}^{\ast }={(X_{im}^T,{\widehat{u}}_i)}^T$, and $f_b$ denotes the density function of the $b_i$’s assumed to be known up to a vector of parameters $\gamma$. Define $L_{im}=max\{U_{ij}:U_{ij}<T_{im},j=0,\dots ,K_i^{\ast }\}$ and $R_{im}=min\{U_{ij}:U_{ij}\ge T_{im},j=1,\dots ,K_i^{\ast }+1\}$, where $U_{i0}=0$ and $U_{i,K_i^{\ast }+1}=\infty$. Then, $(L_{im},R_{im}]$ represents the shortest time interval that brackets $T_{im}$ and the likelihood function above can be rewritten as

$$\begin{array}{ccc}\hfill L_n(\varphi ,\Lambda ,\gamma |u_i^{\prime }s)=& & \prod _{i=1}^n\int \prod _{m=1}^M\left\{exp(-G_m[{\int }_0^{L_{im}}exp(x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i)d{\Lambda }_m\left(t\right)])\right.\hfill \\ & & \left.-exp(-G_m[{\int }_0^{R_{im}}exp(x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i)d{\Lambda }_m\left(t\right)])\right\}{f}_b(b_i\mid \gamma )d{b}_i.\hfill \end{array}$$

By following Huang and Wang (2004) [27] and others, it is natural to estimate $\varphi$ and $\Lambda$ by their values that maximize the approximated likelihood function $L_n(\varphi ,\Lambda ,\gamma |{\widehat{u}}_i^{\prime }s)$.

For the maximization of $L_n(\varphi ,\Lambda ,\gamma |{\widehat{u}}_i^{\prime }s)$, note that it involves the unknown functions ${\Lambda }_m$’s and integrations. To deal with them, for the former, we propose to adopt the nonparametric approach. More specifically, for each $m=1,\dots ,M$, let $0=t_{m0}<t_{m1}<\dots <t_{m{k}_m}<\infty$ denote the ordered sequence of all $L_{im}$ and $R_{im}$ with $R_{im}<\infty$ and assume that ${\Lambda }_m$ is a step function that jumps only at the $t_{mq}$’s with the jump sizes ${\lambda }_{mq}$’s. Then, $L_n(\varphi ,\Lambda ,\gamma |u_i^{\prime }s)$ can be expressed as

$$L_n(\varphi ,\Lambda ,\gamma |u_i^{\prime }s)=\prod _{i=1}^n\int \prod _{m=1}^M\left\{exp\left(-G_m\left[\sum _{t_{mq}\le L_{im}}exp\{x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i\}{\lambda }_{mq}\right]\right)\right.$$

$$\left.-exp\left(-G_m\left[\sum _{t_{mq}\le R_{im}}exp\{x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i\}{\lambda }_{mq}\right]\right)\right\}{f}_b(b_i\mid \gamma )d{b}_i.$$

(3)

In the following, we will develop an EM algorithm for the maximization with the focus on the situation where $f_b$ is a multivariate normal distribution with the covariance matrix $\Sigma \left(\gamma \right)$ depending on the q-dimensional unknown parameter $\gamma$. The algorithm is valid for other distributions and some comments on this will be given below. It is worth to point out that as mentioned above, the idea discussed above has been used by Huang and Wang (2004) [27] and Wang et al. (2016) [23], among others. However, the problem discussed here is different or much more general than the existing literature.

3.2. EM Algorithm

In this subsection, we will develop an EM algorithm for the maximization of $L_n(\varphi ,\Lambda ,$$\gamma \mid {\widehat{u}}_i^{\prime }s)$, and for this, we will first discuss the data augmentation. Let the ${\xi }_{im}$’s denote the random sample of size n from the density $f_m\left(t\right)$. Then, we can rewrite the observed likelihood function as

$$\begin{array}{cc}\hfill L_n& (\varphi ,\Lambda ,\gamma |u_i^{\prime }s)=\prod _{i=1}^n\int \left\{\prod _{m=1}^M{\int }_{{\xi }_{im}}exp\{-{\xi }_{im}\sum _{t_{mq}\le L_{im}}{\lambda }_{mq}exp(x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i)\}\right.\hfill \\ \hfill & \left.{\left[1-exp\{-{\xi }_{im}\sum _{L_{im}<t_{mq}\le R_{im}}{\lambda }_{mq}exp(x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i)\}\right]}^{I(R_{im}<\infty )}{f}_m\left({\xi }_{im}\right)d{\xi }_{im}\right\}\times f_b(b_i\mid \gamma )d{b}_i.\hfill \end{array}$$

(4)

Moreover, let the $W_{imq}$’s denote the random sample of size n from the Poisson distributions with means ${\xi }_{im}{\lambda }_{mq}exp(x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i)$ given ${\xi }_{im}$, and define $A_{im}={\sum }_{t_{mq}\le L_{im}}{W}_{imq}$ and $B_{im}=I(R_{im}<\infty ){\sum }_{L_{im}<t_{mq}\le R_{im}}{W}_{imq}$ such that

$$\begin{array}{cc}\hfill P(A_{im}=0,& B_{im}>0|L_{im},R_{im},x_{im}^{\ast })=exp\{-{\xi }_{im}\sum _{t_{mq}\le L_{im}}{\lambda }_{mq}exp(x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i)\}\hfill \\ \hfill & {\left[1-exp\{-{\xi }_{im}\sum _{L_{im}<t_{mq}\le R_{im}}{\lambda }_{mq}exp(x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i)\}\right]}^{I(R_{im}<\infty )}.\hfill \end{array}$$

It is easy to see that the maximization of (4) is equivalent to maximizing the likelihood function based on the data ($L_{im},R_{im},x_{im}^{\ast },A_{im}=0,B_{im}>0$) $(i=1,\dots ,n;m=1,\dots ,M)$. Based on this, for the development of the EM algorithm, it is natural to use the $W_{imq}$’s, ${\xi }_{im}$’s and $b_i$’s to augment the observed data. As a consequence, one can derive the resulting pseudo complete data log-likelihood function as

$$\begin{array}{cc}\hfill l_c(\varphi ,\Lambda ,\gamma |u_i^{\prime }s)=& \sum _{i=1}^n\left\{\sum _{m=1}^M\left(\sum _{q=1}^{k_m}I(t_{mq}\le R_{im}^{\ast })\left[W_{imq}log\{{\xi }_{im}{\lambda }_{mq}exp(x_{im}^{\ast T}\beta +Z_{im}^T{b}_i)\}\right.\right.\right.\hfill \\ \hfill & \left.\left.-{\xi }_{im}{\lambda }_{mq}exp(x_{im}^{\ast T}\beta +Z_{im}^T{b}_i)-log(W_{imq}!)\right]+log{f}_m\left({\xi }_{im}\right)\right)\hfill \\ \hfill & \left.-\frac{d_i}{2}log\left(2\pi \right)-\frac{1}{2}log\left|{\Sigma }_i\left(\gamma \right)\right|-\frac{b_i^T{\Sigma }_i{\left(\gamma \right)}^{-1}{b}_i}{2}\right\},\hfill \end{array}$$

(5)

where $R_{im}^{\ast }=L_{im}I(R_{im}=\infty )+R_{im}I(R_{im}<\infty )$.

Now, we consider the E-step of the EM algorithm. At the $(s+1)$th iteration and given ${({\varphi }^s,{\Lambda }^s,{\gamma }^s)}^T$, we need to determine

$$Q(\varphi ,\Lambda ,\gamma |{\varphi }^s,{\Lambda }^s,{\gamma }^s)=E\left[l_c(\varphi ,\Lambda ,\gamma |u_i^{\prime }s,\mathcal{O},{\varphi }^s,{\Lambda }^s,{\gamma }^s)\right]$$

$$=\sum _{i=1}^n\left\{\sum _{m=1}^M\left(\sum _{q=1}^{k_m}I(t_{mq}\le R_{im}^{\ast })\left[E\left[W_{imq}log\left\{{\xi }_{im}{\lambda }_{mq}exp(x_{im}^{\ast T}\beta +Z_{im}^T{b}_i)\right\}\right.\right.\right.\right.$$

$$\left.\left.\left.-{\xi }_{im}{\lambda }_{mq}exp(x_{im}^{\ast T}\beta +Z_{im}^T{b}_i)\right]-E\left[log(W_{imq}!)\right]\right]+log{f}_m\left({\xi }_{im}\right)\right)$$

$$\left.-\frac{d_i}{2}log\left(2\pi \right)-\frac{1}{2}log\left|{\Sigma }_i\left(\gamma \right)\right|-\frac{E\left[b_i^T{\Sigma }_i{\left(\gamma \right)}^{-1}{b}_i\right]}{2}\right\}$$

under the multivariate normal distribution with the covariance matrix ${\Sigma }_i\left(\gamma \right)$. To calculate the conditional expectations $E[{\xi }_{im}exp(x_{im}^{\ast T}\beta +Z_{im}^T{b}_i)]$, $E\left[W_{imq}\right]$ and $E\left[b_i^T{\Sigma }^{-1}\left(\gamma \right)b_i\right]$ given the observed data, we need to employ the joint density of ${\xi }_{im}$ and $b_i$ given the observed data, which is proportional to

$$\begin{array}{cc}\hfill & \prod _{m=1}^M\left[exp\left\{-{\xi }_{im}\sum _{t_{mq}\le L_{im}}exp(x_{im}^{\ast T}\beta +Z_{im}^T{b}_i){\lambda }_{mq}\right\}\right.\hfill \\ \hfill & -\left.I(R_{im}<\infty )exp\left\{-{\xi }_{im}\sum _{t_{mq}\le R_{im}}exp(x_{im}^{\ast T}\beta +Z_{im}^T{b}_i){\lambda }_{mq}\right\}\right]\hfill \\ \hfill & \times f_m\left({\xi }_{im}\right){\left(2\pi \right)}^{-d_i/2}{\left|{\Sigma }_i\left(\gamma \right)\right|}^{-1/2}exp\{-\frac{b_i^T{\Sigma }_i{\left(\gamma \right)}^{-1}{b}_i}{2}\}.\hfill \end{array}$$

Note that the conditional expectation of $W_{imq}$ for $t_{mq\le R_{im}^{\ast }}$ given ${\xi }_{im}(m=1,\dots ,M)$, $b_i$ and the observed data is given by

$$\widehat{E}\left(W_{imq}\right|{\xi }_{im},b_i)=I(L_{im}<t_{mq}\le R_{im})\frac{{\lambda }_{mq}{\xi }_{im}exp(x_{im}^{\ast T}\beta +Z_{im}^T{b}_i)}{1-exp\{-{\sum }_{L_{im}<t_{m{q}^{\prime }}\le R_{im}}{\lambda }_{m{q}^{\prime }}{\xi }_{im{q}^{\prime }}exp(x_{im}^{\ast T}\beta +Z_{im}^T{b}_i)\}}.$$

In the M-step, we can employ the Newton–Raphson method to update ${\beta }_m$ based on the equation

$$\sum _{i=1}^n\sum _{m=1}^M\sum _{q=1}^{k_m}I(t_{mq}\le R_{im}^{\ast })\widehat{E}\left(W_{imq}\right)\left\{x_{im}^{\ast }-\frac{{\sum }_{i^{\prime }=1}^{n}I(t_{mq}\le R_{i^{\prime }m}^{\ast })\widehat{E}\{{\xi }_{i^{\prime }m}exp(x_{i^{\prime }m}^{\ast T}{\beta }_m+Z_{i^{\prime }m}^T{b}_{i^{\prime }})\}{x}_{i^{\prime }m}^{\ast }}{{\sum }_{i^{\prime }=1}^{n}I(t_{mq}\le R_{i^{\prime }m}^{\ast })\widehat{E}\{{\xi }_{i^{\prime }m}exp(x_{i^{\prime }m}^{\ast T}{\beta }_m+Z_{i^{\prime }m}^T{b}_{i^{\prime }})\}}\right\}=0.$$

For estimation of ${\lambda }_{mq}$, we have the closed form expression

$${\lambda }_{mq}=\frac{{\sum }_{i=1}^{n}I(t_{mq}\le R_{im}^{\ast })\widehat{E}\left(W_{imq}\right)}{{\sum }_{i=1}^{n}I(t_{mq}\le R_{im}^{\ast })\widehat{E}\{{\xi }_{im}exp(x_{im}^{\ast T}{\beta }_m+Z_{im}^T{b}_i)\}},$$

(6)

for $q=1,\dots ,k_m$ and $m=1,\dots ,M$. To estimate $\gamma$, one can maximize $-log\parallel {\Sigma }_i\left(\gamma \right)\parallel -\widehat{E}\left\{b_i^T{\Sigma }_i^{-1}\left(\gamma \right)b_i\right\}$ with the ${\Sigma }_i$’s given by $\Sigma =n^{-1}{\sum }_{i=1}^n\widehat{E}\left(b_i^T{b}_i\right)$.

Now, we summarize the EM algorithm described above as follows.

Step 1. Choose initial estimates ${\varphi }^{\left(0\right)},{\Lambda }^{\left(0\right)},{\gamma }^{\left(0\right)}$ of $\varphi ,\Lambda ,\gamma$, respectively.

Step 2. In the $(s+1)$th iteration, calculate $\widehat{E}[{\xi }_{im}exp(x_{im}^{\ast T}{\beta }^{\left(s\right)}+Z_{im}^T{b}_i)]$, $\widehat{E}\left[W_{imq}\right]$ and $\widehat{E}[b_i^T{\Sigma }^{-1}$ $\left({\gamma }^{\left(s\right)}\right)b_i]$ by using, for example, the Gaussian quadrature method. $s=0,1,2,\dots$

Step 3. Update ${\Lambda }^{(s+1)}$ by (6) with current ${\varphi }^{\left(s\right)},{\gamma }^{\left(s\right)}$, and then update ${\gamma }^{(s+1)}$ by maximizing $-log\parallel {\Sigma }_i\left(\gamma \right)\parallel -\widehat{E}\left\{b_i^T{\Sigma }_i^{-1}\left(\gamma \right)b_i\right\}$. In addition, estimate ${\varphi }^{(s+1)}$ by employing the one-step Newton–Raphson method.

Step 4. Repeat Steps 2–3 until the convergence such that the absolute difference of the log-likelihood values between two consecutive iterations is less than a given positive value $\epsilon$ such as ${10}^{-3}$.

4. Asymptotic Properties

Let $\widehat{\theta }={\left({\widehat{\beta }}_1^{\mathrm{T}},\dots ,{\widehat{\beta }}_M^{\mathrm{T}},{\widehat{\gamma }}^{\mathrm{T}},{\widehat{\Lambda }}_1,\dots ,{\widehat{\Lambda }}_M\right)}^{\mathrm{T}}$ denote the estimator of $\theta =({\beta }_1^{\mathrm{T}},\dots ,{\beta }_M^{\mathrm{T}},{\gamma }^{\mathrm{T}},$${\Lambda }_1,\dots ,{\Lambda }_M{)}^{\mathrm{T}}$ defined above and ${\theta }_0={({\beta }_{01}^{\mathrm{T}},\dots ,{\beta }_{0M}^{\mathrm{T}},{\gamma }_0^{\mathrm{T}},{\Lambda }_{01},\dots ,{\Lambda }_{0M})}^{\mathrm{T}}$ the true value of $\theta$. Define ${\zeta }_0={\left({\beta }_{01}^{\mathrm{T}},\dots ,{\beta }_{0M}^{\mathrm{T}},{\gamma }_0^{\mathrm{T}}\right)}^{\mathrm{T}}$ and $\widehat{\zeta }={\left({\widehat{\beta }}_1^{\mathrm{T}},\dots ,{\widehat{\beta }}_M^{\mathrm{T}},{\widehat{\gamma }}^{\mathrm{T}}\right)}^{\mathrm{T}}$. In this section, we will establish the asymptotic properties of $\widehat{\theta }$, and for this, we first describe the regularity conditions needed.

Define

$$Q_m^{\ast }\left(t,b;{\beta }_m,{\Lambda }_m\right)=exp\left(-G_m\left[{\int }_0^{t}exp\{{\beta }_m^T{x}_m^{\ast }+b^T{Z}_m\}\right]\right)d{\Lambda }_m\left(s\right),$$

$D_m(U_m,b;{\beta }_m,{\Lambda }_m)={\sum }_{l=0}^{K_m}{\Delta }_{ml}\{Q_m^{\ast }(U_{ml},b;{\beta }_m,{\Lambda }_m)-Q_m^{\ast }(U_{m,l+1},b;{\beta }_m,{\Lambda }_m)\}$, $U_m=$$(U_{m1}$, $\dots ,U_{m,K_m})$, ${\Delta }_{ml}=I(U_{ml}⩽T_m<U_{m,l+1})$, and $p\left(b\right|\gamma )={\left(2\pi \right)}^{-d/2}{|\Sigma \left(\gamma \right)|}^{-1/2}$$exp(-b^{\mathrm{T}}\Sigma {\left(\gamma \right)}^{-1}b/2)$. For the asymptotic properties of $\widehat{\theta }$, we need the following regularity conditions.

Condition 1.

The true value ${\zeta }_0$ belongs to a known compact set $\mathbb{A}⨂\mathbb{B}⨂\mathbb{C}$, where $\mathbb{A}$ denotes a compact set of ${\mathbb{R}}^{\mathrm{pM}}$, $\mathbb{B}$ a compact set in ${\mathbb{R}}^{\mathrm{M}}$, and $\mathbb{C}$ a compact set of ${\mathbb{R}}^{\mathrm{q}}$ in the domain of γ such that $\Sigma \left(\gamma \right)$ is a positive-definite matrix with eigenvalues bounded away from zero and ∞. In addition, the true value ${\Lambda }_{0m}(·)$ is continuously differentiable with positive derivatives in $[0,{\tau }_m]$.

Condition 2.

The covariate vector $X_m$ and $Z_k$ are bounded in $[0,{\tau }_m]$.

Condition 3.

For the transformation function $G_m$, assume that it is twice continuously differentiable on $[0,\infty )$ with $G_m\left(0\right)=0$, $G_m^{{}^{\prime }}\left(x\right)>0$ and $G_m(\infty )=\infty$.

Condition 4.

Assume that ${sup}_{\gamma \in \mathbb{C}}{\int }_{b}g\left(b\right)p^{\left(j\right)}\left(b\right|\gamma )db<\infty$ for any smooth function $g(·)$ and $j=0,1,2$. Here, $p^{\left(j\right)}\left(b\right|\gamma )$ denotes the jth derivative of $p\left(b\right|\gamma )$ with respect to γ.

Condition 5.

If there exists a vector u and some constants $v_m,m=1,\dots ,M$ such that

$$\begin{array}{cc}\hfill & {\left.\left(u^T\frac{\partial }{\partial \zeta }+\sum _{m=1}^M{v}_m\frac{\partial }{\partial y_m}\right)\right|}_{\left(\zeta ,y_1,\dots ,y_M\right)=\left({\zeta }_0,{\Lambda }_{10}\left(c_1\right),\dots ,{\Lambda }_{M0}\left(c_M\right)\right)}\hfill \\ \hfill & ·log{\int }_b\prod _{m=1}^M{D}_m(U_m,b;\beta ,{\beta }_u,{\Lambda }_m)p(b\mid \gamma )\mathrm{d}b=0\hfill \end{array}$$

for each of these values, then $u=0_{pM+M+q}$ and $v_m=0$. In addition, $0_{pM+M+q}$ denotes a $(pM+M+q)$-dimensional vector of zeros.

Condition 6.

Assume that $P({\tau }_m\ge {\tau }_0,exp\left(u\right)>0)>0$ for the follow-up time ${\tau }_m$ and latent variable u, where ${\tau }_0$ denotes the longest study time and the variance of $exp\left(u\right)$ is bounded and there exists a positive small constant $ϵ>0$ such that $exp\left(u\right)>ϵ$ almost surely. Moreover, for ${\tau }_m$ and u, the function $F\left(s\right)=E[exp\left(u\right)I({\tau }_m\ge s)]$ is continuous for $s\in [0,{\tau }_0]$.

Note that Conditions 1 and 2 are standard conditions in survival analysis, and it is easy to check that Condition 3 on the transformation function holds for the logarithmic family $G_r\left(x\right)=r^{-1}log(1+rx)(r⩾0)$ and the Box–Cox family $G_d\left(x\right)=d^{-1}\left\{{(1+x)}^d-1\right\}(d⩾0)$. Moreover, Condition 4 holds for modeling multivariate data with frailty models, and Condition 5 is required for the identifiability of the model. In addition, Condition 6 describes the relationship between the latent variable u and the parameters of interest. Most of the conditions above are purely for technique purposes and hold in general, in particular, for periodic follow-up studies.

Let $\parallel ·\parallel$ denote the Euclidean norm and define $\mathbb{P}f=\int f\left(x\right)\mathrm{d}P\left(x\right)$ and ${\mathbb{P}}_{n}f=n^{-1}{\sum }_{i=1}^{n}f\left(X_i\right)$ for a function f and a random variable X with distribution P. The following two theorems give the asymptotic properties of $\widehat{\theta }$.

Theorem 1.

Suppose that Conditions 1–6 hold. Then, as $n\to \infty$, we have that $∥\widehat{\zeta }-{\zeta }_0∥$$+{\sum }_{m=1}^M{sup}_{t\in \left[0,{\tau }_m\right]}\left|{\widehat{\Lambda }}_m\left(t\right)-{\Lambda }_{0m}\left(t\right)\right|\to 0$ almost surely.

Theorem 2.

Suppose that Conditions 1–6 hold. Then, as $n\to \infty$, we have that $\sqrt{n}({\widehat{\zeta }}_n-{\zeta }_0){\to }^{d}N(0,I_0^{-1})$, where $I_0=\mathbb{P}\left\{\tilde{l}\left({\theta }_0\right)\tilde{l}{\left({\theta }_0\right)}^T\right\}$ with $\tilde{l}\left({\theta }_0\right)$ given in the Appendix A.

We will sketch the proof for the results described above in Appendix A. For inference about $\zeta$, it is apparent that one needs to estimate the covariance matrix, and for this, one can see from Appendix A that it would be difficult to derive a consistent estimator of $I_0$. Thus, we propose to employ the profile likelihood approach to estimate the covariance matrix of $\widehat{\zeta }$ (Murphy & van der Vaart, 2000) [28]. Specifically, let $\mathcal{C}$ denote the set of all step functions with nonnegative jumps at $t_{mq}$ and define ${\mathrm{pl}}_n\left(\zeta \right)={max}_{\Lambda \in \mathcal{C}}log{L}_n(\zeta ,\Lambda )$, the profile log-likelihood. Then, one can estimate the covariance matrix of $\widehat{\zeta }$ by the negative inverse of the matrix with the $(j,k)$th element given by

$$\frac{{\mathrm{pl}}_n\left(\widehat{\zeta }\right)-{\mathrm{pl}}_n\left(\widehat{\zeta }+h_n{e}_k\right)-{\mathrm{pl}}_n\left(\widehat{\zeta }+h_n{e}_j\right)+{\mathrm{pl}}_n\left(\widehat{\zeta }+h_n{e}_k+h_n{e}_j\right)}{h_n^2}.$$

In the above, $e_j$ denotes the jth canonical vector in ${\mathcal{R}}^d$ and $h_n$ is a constant of order $n^{-1/2}$. Note that to calculate $p{l}_n\left(\zeta \right)$ for each $\zeta$, one can reuse the proposed EM algorithm with $\beta$ held fixed and the only step in the EM algorithm is to explicitly evaluate $\widehat{E}\left(W_{imq}\right)$ and $\widehat{E}\left({\xi }_{im}\right)$ to update ${\lambda }_m$ using above. The iteration converges quickly in general by setting ${\widehat{\lambda }}_m$ to be the initial value.

5. A Simulation Study

In this section, we give some of the simulation results obtained from a study performed to evaluate the finite sample performance of the proposed method with the focus on estimation of the ${\beta }_m$’s. In the study, we considered the situation with $M=2$ correlated failure times of interest and two covariates. For the covariates, it was assumed that the first covariate follows the Bernoulli distribution with the success probability of 0.5 and the second covariate, the uniform distribution over $(0,1)$. To generate the true failure times, we first set $Z_{im}$ to be one and generated the latent variables $b_i$’s from the normal distribution $N(0,{\sigma }^2)$ with ${\sigma }^2=0.25$ and the latent variables $u_i$’s from the normal distribution with the mean 0 and variance 1. Then, given the $X_{im}$’s, $Z_{im}$’s, $b_i$’s and $u_i$’s, the $T_{i1}$’s and $T_{i2}$’s were generated under model (1) with $G_m\left(x\right)=r_m^{-1}log(1+r_{m}x)$, ${\Lambda }_1\left(t\right)=log(1+0.5t)$ and ${\Lambda }_2\left(t\right)=0.5t$ for $r_1=r_2=0$, $r_1=r_2=0.5$ or $r_1=r_2=1$, respectively.

For the generation of the observation process and the observed data, we first assumed that the ${\tau }_i^{\prime }s$ follow the uniform distribution over the interval [2, 3] and generated the $K_i^{\ast }$’s from the Poisson distribution with the mean

$${\Lambda }_{ih}({\tau }_i;X_i,u_i)={\tau }_{i}exp(X_i^T\alpha +u_i)$$

given the $X_i$’s and $u_i$’s. Note that in the above, we took ${\Lambda }_{0h}\left(t\right)=t$ and ${\alpha }_m=1$. Given the $K_i^{\ast }$’s, we took $U_{i1}<\dots <U_{i{K}_i^{\ast }}$ to be the order statistics of the random sample of size $K_i^{\ast }$ from the uniform distribution over $(0,{\tau }_i)$. In the following, we considered two sets of true values, ${(0,0,0)}^T$ and ${(0.5,0.5,0.5)}^T$, for the regression parameters ${\beta }_1={({\beta }_{x11},{\beta }_{x12},{\beta }_{u1})}^T$ and ${\beta }_2={({\beta }_{x21},{\beta }_{x22},{\beta }_{u2})}^T$, corresponding to $T_{i1}$ and $T_{i2}$, respectively. The results given below are based on $n=200$ or 400 with 1000 replications.

Table 1. Simulation results on estimation of $\beta$ with the $b_i$’s generated from the normal distribution.

	True Value	${\mathit{r}}_1={\mathit{r}}_2=0$				${\mathit{r}}_1={\mathit{r}}_2=0.5$				${\mathit{r}}_1={\mathit{r}}_2=1$
		Bias	SSE	SEE	CP	Bias	SSE	SEE	CP	Bias	SSE	SEE	CP
$n=200$
${\beta }_{x11}$	0	0.030	0.248	0.252	0.945	0.006	0.289	0.288	0.951	0.031	0.251	0.259	0.945
${\beta }_{x12}$	0	0.028	0.412	0.422	0.957	0.021	0.479	0.476	0.953	0.029	0.418	0.452	0.957
${\beta }_{u1}$	0	−0.024	0.141	0.143	0.938	−0.026	0.164	0.163	0.945	−0.024	0.142	0.142	0.938
${\beta }_{x21}$	0	0.044	0.236	0.236	0.953	0.001	0.278	0.276	0.948	0.044	0.242	0.252	0.953
${\beta }_{x22}$	0	0.001	0.390	0.422	0.949	0.011	0.459	0.460	0.949	0.001	0.403	0.429	0.949
${\beta }_{u2}$	0	−0.013	0.138	0.135	0.960	−0.026	0.159	0.158	0.945	−0.014	0.138	0.140	0.960
${\beta }_{x11}$	0.5	0.040	0.233	0.246	0.944	0.024	0.291	0.290	0.955	0.033	0.314	0.317	0.950
${\beta }_{x12}$	0.5	0.025	0.378	0.410	0.955	−0.007	0.485	0.501	0.951	0.020	0.508	0.533	0.948
${\beta }_{u1}$	0.5	0.010	0.135	0.142	0.956	−0.026	0.165	0.168	0.942	−0.022	0.177	0.179	0.948
${\beta }_{x21}$	0.5	0.028	0.227	0.238	0.949	0.017	0.280	0.282	0.949	0.062	0.307	0.316	0.947
${\beta }_{x22}$	0.5	0.011	0.358	0.436	0.953	0.011	0.466	0.484	0.956	0.054	0.489	0.498	0.949
${\beta }_{u2}$	0.5	−0.005	0.135	0.131	0.943	−0.027	0.160	0.160	0.950	−0.022	0.176	0.179	0.958
$n=400$
${\beta }_{x11}$	0	0.009	0.171	0.178	0.952	0.010	0.200	0.204	0.954	−0.001	0.226	0.228	0.949
${\beta }_{x12}$	0	0.001	0.280	0.283	0.952	−0.003	0.327	0.329	0.954	0.017	0.370	0.374	0.946
${\beta }_{u1}$	0	−0.023	0.097	0.104	0.949	−0.036	0.112	0.113	0.944	0.036	0.127	0.130	0.942
${\beta }_{x21}$	0	0.009	0.162	0.180	0.941	0.007	0.191	0.194	0.946	0.009	0.219	0.222	0.942
${\beta }_{x22}$	0	0.001	0.260	0.287	0.948	0.014	0.314	0.321	0.951	0.003	0.358	0.366	0.955
${\beta }_{u2}$	0	−0.025	0.093	0.096	0.942	−0.037	0.108	0.109	0.948	−0.036	0.124	0.123	0.946
${\beta }_{x11}$	0.5	0.005	0.158	0.163	0.949	0.012	0.187	0.190	0.949	0.031	0.215	0.22	0.942
${\beta }_{x12}$	0.5	0.019	0.253	0.254	0.954	0.023	0.299	0.305	0.957	0.010	0.344	0.363	0.949
${\beta }_{u1}$	0.5	−0.025	0.090	0.096	0.935	−0.031	0.105	0.105	0.937	−0.034	0.121	0.120	0.944
${\beta }_{x21}$	0.5	0.012	0.153	0.166	0.952	0.008	0.183	0.185	0.942	0.024	0.211	0.219	0.946
${\beta }_{x22}$	0.5	0.035	0.239	0.259	0.956	0.007	0.289	0.302	0.949	0.011	0.333	0.348	0.953
${\beta }_{u2}$	0.5	−0.024	0.089	0.094	0.941	−0.032	0.104	0.104	0.94	−0.03	0.120	0.115	0.943

gives the results on the estimation of ${\beta }_1$ and ${\beta }_2$ given by the proposed estimation procedure with $r_1=r_2=0$, $r_1=r_2=0.5$ and $r_1=r_2=1$. Here, we calculated the estimated bias (Bias) given by the average of the estimates minus the true value, the sample standard error (SSE) of the estimates, the average of the estimated standard errors (ESE) and the $95\%$ empirical coverage probability (CP). The results suggest that the proposed estimator of the regression parameters seems to be unbiased and the variance estimation based on the profile likelihood approach also seems to be reasonable. Furthermore, the results on the empirical coverage probabilities indicate that the normal approximation to the distribution of the proposed estimator of the regression parameters appears to be appropriate. In addition, the results got better in general with the increasing sample size, as expected.

As mentioned before, the proposed estimation procedure can be applied to any distribution for the latent variables $b_i$’s. To see this, we repeated the study above, except that we generated the $b_i$’s from the uniform distribution over $(-1,1)$, and

Table 2. Simulation results on estimation of $\beta$ with the $b_i$’s generated from the uniform distribution and $r_1=r_2=0$.

	True Value	Bias	SSE	ESE	CP
${\beta }_{x11}$	0	−0.005	0.238	0.240	0.948
${\beta }_{x12}$	0	−0.004	0.404	0.404	0.953
${\beta }_{u1}$	0	0.001	0.134	0.134	0.946
${\beta }_{x21}$	0	−0.001	0.223	0.232	0.950
${\beta }_{x22}$	0	−0.004	0.374	0.412	0.945
${\beta }_{u2}$	0	0.003	0.126	0.126	0.957
${\beta }_{x11}$	0.5	0.031	0.223	0.231	0.948
${\beta }_{x12}$	0.5	0.035	0.361	0.401	0.956
${\beta }_{u1}$	0.5	−0.021	0.129	0.134	0.947
${\beta }_{x21}$	0.5	0.022	0.215	0.231	0.948
${\beta }_{x22}$	0.5	0.054	0.342	0.412	0.947
${\beta }_{u2}$	0.5	−0.012	0.129	0.129	0.955

presents the obtained results on the estimation of ${\beta }_1$ and ${\beta }_2$ with $n=200$ and $r_1=r_2=0$. One can see that they are similar to those given in Table 1 and again suggest that the proposed approach seems to work well for the situations considered. To see the performance of the proposed approach with different types of covariates, we also repeated the study giving Table 1, except that both covariates were assumed to follow the standard normal distribution and give the obtained results with $n=200$ and $r_1=r_2=0$ in

Table 3. Simulation results on estimation of $\beta$ with the covariates generated from the normal distribution and $r_1=r_2=0$.

	True Value	Bias	SSE	ESE	CP
${\beta }_{x11}$	0	−0.003	0.127	0.125	0.949
${\beta }_{x12}$	0	0.003	0.126	0.126	0.953
${\beta }_{u1}$	0	−0.006	0.142	0.145	0.950
${\beta }_{x21}$	0	0.003	0.120	0.119	0.938
${\beta }_{x22}$	0	−0.001	0.120	0.119	0.941
${\beta }_{u2}$	0	0.004	0.133	0.141	0.957
${\beta }_{x11}$	0.5	0.037	0.142	0.141	0.946
${\beta }_{x12}$	0.5	0.027	0.141	0.141	0.954
${\beta }_{u1}$	0.5	−0.011	0.146	0.152	0.948
${\beta }_{x21}$	0.5	0.037	0.137	0.135	0.945
${\beta }_{x22}$	0.5	0.033	0.138	0.134	0.939
${\beta }_{u2}$	0.5	−0.008	0.143	0.148	0.949

. They indicate that the proposed estimation procedure seems to be robust to different types of covariates.

Note that in the proposed estimation procedure, it has been assumed that the observation process ${\tilde{N}}_i\left(t\right)$ is a non-homogeneous Poisson process and one may be interested in the performance of the proposed method if this assumption is not true. To see this, we repeated the study giving Table 1, except that the ${\tilde{N}}_i\left(t\right)$’s were assumed to be mixed Poisson processes with the intensity function

$${\lambda }_{ih}\left(t\right|X_i,u_i)=v_i{\lambda }_{0h}\left(t\right)exp\left(X_i^T\alpha +u_i\right)$$

given the $v_i$’s, where the $v_i$’s were generated from the gamma distribution.

Table 4. Simulation results on estimation of $\beta$ with mixed Poisson observation processes and $r_1=r_2=0$.

	True Value	Bias	SSE	ESE	CP
${\beta }_{x11}$	0	0.002	0.246	0.250	0.952
${\beta }_{x12}$	0	0.013	0.418	0.429	0.951
${\beta }_{u1}$	0	0.001	0.137	0.136	0.951
${\beta }_{x21}$	0	−0.002	0.233	0.241	0.948
${\beta }_{x22}$	0	0.006	0.387	0.419	0.946
${\beta }_{u2}$	0	−0.001	0.130	0.131	0.949
${\beta }_{x11}$	0.5	0.037	0.233	0.241	0.948
${\beta }_{x12}$	0.5	0.075	0.376	0.400	0.949
${\beta }_{u1}$	0.5	−0.034	0.133	0.136	0.945
${\beta }_{x21}$	0.5	0.037	0.226	0.257	0.953
${\beta }_{x22}$	0.5	0.038	0.357	0.435	0.949
${\beta }_{u2}$	0.5	−0.029	0.132	0.142	0.952

presents the results on the estimation of ${\beta }_1$ and ${\beta }_2$ given by the proposed approach with $n=200$ and $r_1=r_2=0$, and they indicate that the approach seems to be robust with the processes ${\tilde{N}}_i\left(t\right)$’s.

For the initial value in the EM algorithm here, we set ${\varphi }^{\left(0\right)}=0$, ${\lambda }_{mq}^{\left(0\right)}=\frac{1}{k_m}$, $q=1,\dots ,k_m$, and ${\gamma }^{\left(0\right)}=0.25$. It is worth to point out that we did try other initial values and the proposed EM algorithm seems to be robust with respect to the selection of the initial values. In other words, we did not encounter non-convergence issue in the simulation study. We also considered some other setups, including multivariate cases and the case with more than one covariate and obtained similar results.

6. An Application

In the section, the estimation procedure proposed in the previous sections is applied to the set of bivariate interval-censored data arising from an AIDS clinical trial, AIDS Clinical Trial Group 181, described in Goggins and Finkelstein (2000) [11]. The study concerns the opportunistic infection cytomegalovirus (CMV) and examined the study patients periodically. At each clinical visit or observation, among other information, the blood and urine samples were collected and tested to detect the existence of the CMV virus in the sample, which is also commonly referred to as the shedding of the virus. In addition, for each patient, the CD4 count, indicating the status of a person’s immune system and being commonly used to measure the stage of HIV infection, was also recorded at the entry time. For the analysis here, we are mainly interested in if the baseline CD4 account, the indicator of the initial stage of HIV disease, is related to the CMV shedding risk in either blood or urine.

The data set consists of 204 subjects, and they belong to two groups based on their baseline CD4 counts, either less than 75 or otherwise. More specifically, the two groups have 111 and 93 patients, respectively. On the observation of the CMV shedding times, some patients gave left-censored observations and some right-censored observations. The others provided some intervals or interval-censored observations, given by the last negative and first positive test dates. That is, we have bivariate interval-censored data on the CMV shedding times in the blood and urine. The percentages of right-censored observations for the CMV shedding times in the blood and urine are about $85\%$ and $43\%$, respectively, which indicate that the CMV shedding risk in the urine may be higher than that in the blood. For the application of the proposed estimation procedure, let $T_{i1}$ and $T_{i2}$ denote the CMV shedding times in the blood and urine associated with the ith patient, respectively, and define $X_i=1$ if the ith subject’s baseline CD4 count was less than 75 and 0 otherwise. As in the simulation study, we took $G_m\left(x\right)=r_m^{-1}log(1+r_{m}x)$ and set $Z_i=1$ for all patients.

Table 5. Analysis results for the AIDS clinical trials data.

${\mathit{r}}_1$	${\mathit{r}}_2$	${\mathit{\beta }}_{\mathit{blood}}$	${\mathit{\beta }}_{\mathit{urine}}$	${\mathit{SE}}_{\mathit{blood}}$	${\mathit{SE}}_{\mathit{urine}}$	${\mathit{P}}_{\mathit{blood}}$	${\mathit{P}}_{\mathit{urine}}$	AIC	BIC
0	0	2.312	1.143	0.396	0.142	0.000	0.001	727.917	744.507
	0.5	2.333	1.291	0.401	0.208	0.000	0.002	733.972	750.563
	1	2.381	1.437	0.413	0.256	0.000	0.002	738.959	755.550
0.5	0	2.483	1.141	0.424	0.141	0.000	0.001	727.327	743.918
	0.5	2.504	1.288	0.430	0.206	0.000	0.002	733.405	749.996
	1	2.552	1.435	0.442	0.255	0.000	0.002	738.426	755.016
1	0	2.652	1.140	0.451	0.140	0.000	0.001	726.826	743.416
	0.5	2.670	1.285	0.457	0.205	0.000	0.002	732.924	749.514
	1	2.717	1.433	0.469	0.254	0.000	0.002	737.981	754.572

presents the estimation results given by the proposed approach for different combinations of $r_1=0,0.5,1$ and $r_2=0,0.5,1$, and they include the estimated covariate effects, ${\widehat{\beta }}_{blood}$ and ${\widehat{\beta }}_{urine}$, the estimated standard errors (SE) and the p-values for testing no covariate effect (P). In addition, we have calculated the Akaike Information Criterion (AIC, Akaike, 1973 [29]) and Bayesian Information Criterion (BIC, Schwarz, 1978 [30]) for the selection of the optimal model. One can see from the table that the AIC and BIC values are quite close for all combinations of $r_1$ and $r_2$, and the same is true for the estimated effects. By choosing $r_1=r_2=0$, which correspond to the proportional hazards models for both the $T_{i1}$’s and $T_{i2}$’s, we have ${\widehat{\beta }}_{blood}=2.312$ and ${\widehat{\beta }}_{urine}=1.143$ with the estimated standard errors being $0.396$ and $0.142$, respectively. They suggest that the patients with lower CD4 at the baseline experienced CMV shedding in both blood and urine significantly early. To provide a graphical view about the difference between the CMV shedding in the blood and urine, Figure 1 presents the estimates of the baseline marginal survival functions given by the proposed method with $r_1=r_2=0$ for the CMV shedding times in the blood and urine, respectively. They suggest that as discussed above, the CMV shedding in the urine occurred much earlier than in the blood.

In addition, with $r_1=r_2=0$, the proposed method yielded ${\widehat{\beta }}_{u1}=2.845$ and ${\widehat{\beta }}_{u2}=0.958$ with the estimated standard errors of $0.112$ and $0.116$, respectively. They indicate that the observation process was significantly correlated with the CMV shedding times in both blood and urine. That is, we had dependent or informative censoring. To investigate the effects of informative censoring on the covariate effects, we assumed that ${\beta }_{u1}={\beta }_{u2}=0$, meaning independent interval censoring, and obtained ${\widehat{\beta }}_{blood}=1.560$ and ${\widehat{\beta }}_{urine}=1.306$ with the estimated standard errors being $0.514$ and $0.326$, respectively. They would correspond to the p-values of $0.015$ and $0.011$ for testing ${\beta }_{blood}=0$ and ${\beta }_{urine}=0$, respectively. Although these results are similar to those given above, it is apparent that they underestimated the effects of the baseline CD4 on the risks of the CMV shedding times.

7. Discussion and Concluding Remarks

In the preceding sections, the regression analysis of case K multivariate interval-censored failure-time data was discussed under a general class of semiparametric transformation models in the presence of informative censoring. For the problem, an approximate maximum likelihood estimation procedure was proposed and the resulting estimators of the regression parameters were shown to be consistent and asymptotically normal. In the method, the frailty approach was employed to characterize the informative censoring as well as the relationship among the correlated failure times of interest. To implement the proposed approach, a novel EM algorithm was developed and the numerical studies indicated that the proposed method works well in practical situations. In addition, it was applied to a set of real bivariate interval-censored data arising from an AIDS clinical trial.

The proposed approach can be seen as a generalization of the method given by Zeng et al. (2017) [15] to allow for informative interval censoring, which can occur quite often, as discussed above and in the literature. In particular, it has been shown that in the presence of informative censoring, the analysis that ignores it could lead to biased or misleading results and conclusions. The proposed method has the advantages that it does not need or impose an assumption on the distribution of the latent variables and it is quite flexible and can be easily implemented. Moreover, the type of the data considered here includes most types of the failure-time data discussed in the literature as special cases and the model (1) gives many commonly used models.

As discussed above, although model (1) is quite flexible, it may not be straightforward to choose an optimal model for a given set of data, and one commonly used procedure for this is to apply the AIC or BIC. As an alternative, one may prefer to develop a model-checking or data-driven technique. However, this may be difficult and such a method does not seem to exist even for simple types of multivariate interval-censored data. It is worth noting that instead of the proposed approximation maximum likelihood estimation method, one may consider a full maximum likelihood estimation procedure. For this, however, one would need to specify or postulate some distributions for the latent variables $b_i$’s, which may be hard to be verified, and also the implementation would be much more complicated.