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

: 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.


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.
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 noninformative 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 intervalcensored 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 intervalcensored 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.

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 < . . .< U iK * i and a followup or stopping time τ i , where K * i denotes the number of potential observations, i = 1, . . ., 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 Ñi (t) = ∑ K * i j=1 I(U ij ≤ t), describing the observation process on the subject that jumps only at the observation times, i = 1, . . ., 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 Ñi (t), 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 Here, G m (.) is a known strictly increasing transformation function, Λ m (.) is an unknown baseline cumulative hazard function, Z im contains 1 and part of the covariates X im , and β m = (β T xm , β um ) T denotes the vector of unknown regression parameters.For the observation process, it will be assumed that Ñi (t) is a non-homogeneous Poisson process satisfying for the intensity function given X im and u i .Here, λ 0h is an unknown continuous baseline intensity function, X i = (X T i1 , . . ., X T iM ) T , and α T = (α T 1 , . . ., α T M ), which is a vector of regression parameters as β m .In the following, it will be assumed that given X im , b i and u i , T i1 , . . ., T iM are independent, and given X im and u i , T i and Ñi are independent.Moreover, τ i is independent of T i and Ñ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 β um denotes the degree of the correlation between the failure-time process and the observation process, and they will be independent if β 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 frailtyinduced transformations In the above, f m (t) denotes the density function of a frailty variable with support [0, ∞].By setting f m (t) to be the gamma density with mean 1 and variance r m , it gives the class of logarithmic transformations G m (x) = r −1 m 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 (x) = log(1 + x) and gives the proportional hazards model with To describe the observed data, define In the following, it will be assumed that the observed data have the form where That is, we observe case K interval-censored data (Sun, 2006).

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 φ T = (β T 1 , . . ., β T M ), 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 follows the Poisson distribution with the mean given X i and u i , i = 1, . . ., n.The other is that the observation times U i1 , . . ., U iK * i can be seen as the order statistics of a set of i.i.d.random variables with the density function One can see that the function above does not depend on neither X i nor u i , which suggests that the function Λ 0h (t) can be estimated by ) .
In the above, the s (l) 's denote the ordered and distinct values of observation times {U ik }, d (l) the number of the observation times equal to s (l) , and R (l) the number of observation times satisfying U ik ≤ s (l) ≤ τ i among all subjects.Under the assumptions above, it is easy to show that and a class of estimating equations for estimation of α m , m = 1, . . ., M with the ω i 's being some weights.Let αm denote the estimator of α m given by the estimating equations above, which suggests that one can naturally estimate u i by ûi = log . Now, consider estimation of φ 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 where Λ = (Λ 1 , . . ., Λ M ), x * im = (X T im , ûi ) T , and f b denotes the density function of the b i 's assumed to be known up to a vector of parameters γ.
Then, (L im , R im ] represents the shortest time interval that brackets T im and the likelihood function above can be rewritten as By following Huang and Wang (2004) [27] and others, it is natural to estimate φ and Λ by their values that maximize the approximated likelihood function L n (φ, Λ, γ | û i s).
For the maximization of L n (φ, Λ, γ | û i s), note that it involves the unknown functions Λ m 's and integrations.To deal with them, for the former, we propose to adopt the nonparametric approach.More specifically, for each m = 1, . . ., M, let 0 = t m0 < t m1 < . . .< t mk m < ∞ denote the ordered sequence of all L im and R im with R im < ∞ and assume that Λ m is a step function that jumps only at the t mq 's with the jump sizes λ mq 's.Then, L n (φ, Λ, γ | u i s) can be expressed as 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 Σ(γ) depending on the q-dimensional unknown parameter γ.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.

EM Algorithm
In this subsection, we will develop an EM algorithm for the maximization of L n (φ, Λ, γ | û i s), and for this, we will first discuss the data augmentation.Let the ξ im 's denote the random sample of size n from the density f m (t).Then, we can rewrite the observed likelihood function as Moreover, let the W imq 's denote the random sample of size n from the Poisson distributions with means .
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 , A im = 0, B im > 0) (i = 1, . . ., n; m = 1, . . ., M).Based on this, for the development of the EM algorithm, it is natural to use the W imq 's, ξ 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 where . Now, we consider the E-step of the EM algorithm.At the (s + 1)th iteration and given (φ s , Λ s , γ s ) T , we need to determine under the multivariate normal distribution with the covariance matrix Σ i (γ).To calculate the conditional expectations given the observed data, we need to employ the joint density of ξ im and b i given the observed data, which is proportional to Note that the conditional expectation of W imq for t mq≤R * im given ξ im (m = 1, . . ., M), b i and the observed data is given by In the M-step, we can employ the Newton-Raphson method to update β m based on the equation For estimation of λ mq , we have the closed form expression for q = 1, . . ., k m and m = 1, . . ., M. To estimate γ, one can maximize ).Now, we summarize the EM algorithm described above as follows.
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 ε such as 10 −3 .

Asymptotic Properties
Let θ = βT 1 , . . ., βT M , γT , Λ1 , . . ., ΛM T denote the estimator of θ = (β T 1 , . . ., β T M , γ T , Λ 1 , . . . ,Λ M ) T defined above and In this section, we will establish the asymptotic properties of θ, and for this, we first describe the regularity conditions needed.Define For the asymptotic properties of θ, we need the following regularity conditions.Condition 1.The true value ζ 0 belongs to a known compact set A B C, where A denotes a compact set of R pM , B a compact set in R M , and C a compact set of R q in the domain of γ such that Σ(γ) is a positive-definite matrix with eigenvalues bounded away from zero and ∞.In addition, the true value Λ 0m (•) is continuously differentiable with positive derivatives in [0, τ m ].
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(τ m ≥ τ 0 , exp(u) > 0) > 0 for the follow-up time τ m and latent variable u, where τ 0 denotes the longest study time and the variance of exp(u) is bounded and there exists a positive small constant > 0 such that exp(u) > almost surely.Moreover, for τ m and u, the function 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 . 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 • denote the Euclidean norm and define P f = f (x)dP(x) and P n f = n −1 ∑ n i=1 f (X i ) for a function f and a random variable X with distribution P. The following two theorems give the asymptotic properties of θ.
We will sketch the proof for the results described above in Appendix A. For inference about ζ, 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 ζ (Murphy & van der Vaart, 2000) [28].Specifically, let C denote the set of all step functions with nonnegative jumps at t mq and define pl n (ζ) = max Λ∈C log L n (ζ, Λ), the profile log-likelihood.Then, one can estimate the covariance matrix of ζ by the negative inverse of the matrix with the (j, k)th element given by pl In the above, e j denotes the jth canonical vector in R d and h n is a constant of order n −1/2 .Note that to calculate pl n (ζ) for each ζ, one can reuse the proposed EM algorithm with β held fixed and the only step in the EM algorithm is to explicitly evaluate Ê(W imq ) and Ê(ξ im ) to update λ m using above.The iteration converges quickly in general by setting λm to be the initial value.

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 β 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, σ 2 ) with σ 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) For the generation of the observation process and the observed data, we first assumed that the τ i s follow the uniform distribution over the interval [2,3] and generated the K * i 's from the Poisson distribution with the mean given the X i 's and u i 's.Note that in the above, we took Λ 0h (t) = t and α m = 1.Given the K * i 's, we took U i1 < . . .< U iK * i to be the order statistics of the random sample of size K * i from the uniform distribution over (0, τ 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 β 1 = (β x11 , β x12 , β u1 ) T and β 2 = (β x21 , β x22 , β 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 gives the results on the estimation of β 1 and β 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 presents the obtained results on the estimation of β 1 and β 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.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 Ñi (t) 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 Ñi (t)'s were assumed to be mixed Poisson processes with the intensity function given the v i 's, where the v i 's were generated from the gamma distribution.Table 4 presents the results on the estimation of β 1 and β 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 Ñi (t)'s.For the initial value in the EM algorithm here, we set φ (0) = 0, λ (0) mq = 1 k m , q = 1, . . ., k m , and γ (0) = 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.

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 (x) = r −1 m log(1 + r m x) and set Z i = 1 for all patients.Table 5 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, βblood and β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 βblood = 2.312 and β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 βu1 = 2.845 and β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 β u1 = β u2 = 0, meaning independent interval censoring, and obtained βblood = 1.560 and β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 β blood = 0 and β 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.

Discussion and Concluding Remarks
In the preceding sections, the regression analysis of case K multivariate intervalcensored 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 modelchecking 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.
Then, the class of function D m (U m , b; β m , Λ m ) is a Donsker class.By condition 4, we know that l(θ) is bounded away from zero.Therefore, l(θ, O) belongs to some Donsker class due to the preservation property of the Donsker class under the Lipschitz-continuous transformations.Then, we can conclude that sup θ∈M |P n l(θ, O) − Pl(θ 0 , O)| converges in probability to 0 as n → ∞.
By following similar calculations in Section 3 of Chang et al. (2007) [32], we can establish the existence of h * m in the above equation.The efficient Fisher information matrix I 0 for ζ at (ζ 0 , Λ 0 ) is defined as P l(ζ 0 , Λ 0 ) l(ζ 0 , Λ 0 ) T .In the following, we will show that I 0 is positive definite.If the I 0 is singular, then there exists a nonzero vector ν ∈ R (pM+M+q)  such that ν T I 0 ν = 0.It follows that, with probability one, the score function along the submodel ζ 0 + ν, Λ 10 + ν T h * 1 , . . ., Λ M0 + ν T h * M is zero.Therefore, Using Condition 5, we know that ν = 0, and this is a contradiction.Therefore, we can conclude that ν T I 0 ν = 0 implies ν = 0.That is, the efficient Fisher information matrix is positive.Define We know from Example 19.11 of Van der Vaart (1998) [31] that the class of uniformly bounded functions with bounded variations is a Donsker class.By using Theorem 2.10.6 of Van der Vaart and Wellner (1996) [33], we can verify that l(ζ, Λ) is a uniformly bounded Donsker class.In addition, we have proved that θn is consistent.Therefore, This completes the proof of Theorem 2.

Figure 1 .
Figure 1.Estimated marginal survival functions for the CMV shedding times in the blood and urine.

Table 1 .
Simulation results on estimation of β with the b i 's generated from the normal distribution.

Table 2 .
Simulation results on estimation of β with the b i 's generated from the uniform distribution and r 1 = r 2 = 0.

Table 3 .
Simulation results on estimation of β with the covariates generated from the normal distribution and r 1 = r 2 = 0.

Table 4 .
Simulation results on estimation of β with mixed Poisson observation processes and r 1 = r 2 = 0.

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