Bayesian and hierarchical Bayesian analysis of response-time data with concomitant variables

This paper considers the Bayes and hierarchical Bayes approaches for analyzing clinical data on response times with available values for one or more concomitant variables. Response times are assumed to follow simple exponential distributions, with a different parameter for each patient. The analyses are carried out in case of progressive censoring assuming squared error loss function and gamma distribution as priors and hyperpriors. The possibilities of using the methodology in more general situations like doseresponse modeling have also been explored. Bayesian estimators derived in this paper are applied to lung cancer data set with concomitant variables.


INTRODUCTION
In biomedical studies, a considerable interest is laid upon developing statistical techniques for analyzing survival data which utilize information available on concomitant variables.In classical analysis of survival data, several models [1][2][3][4][5][6][7] are used for such situations.The usual proportional hazards (PH) regression model proposed by Cox [8] has been extensively discussed in the literature.Byar et al. [9] and Greenberg et al. [10] presented analysis of survival data assuming linear hazard model in classical set-up.
Bhattacharya et al. [11] discussed for the first time the problem on estimation of survival probabilities adjusting the effect of a single concomitant variable in the Bayesian framework.The present paper presents the Bayesian and hierarchical Bayesian analysis of response-time data in more general situations of more than one concomitant variables available for their effects to be adjusted.The exponential survival model representing the death density function (DDF) corresponding to the survival time Y is assumed.We also assume that the hazard λ for a patient under clinical investigation is linearly related to measurements on 'p' concomitant variables x 1 , x 2 , . . ., x p as follows where β 0 , β 1 , . . . ., β p are (p + 1) unknown parameters and x 0 = 1 is a dummy variable which is set equal to 1 for all individuals for notational symmetry.In (2) β 0 can be interpreted as the underlying hazard rate or the intercept.Of course, it is necessary that the right hand side of (2) be positive.The above hazard model can also be written as where x  = (x 0 , x 1 , . . . ., x p ) is a 1 × (p + 1) vector of concomitant variables measured on the individual under clinical investigation and A natural extension of the model ( 2) is a dose-response model with hazard as a polynomial function of concomitant variables covering the situations wherein some concomitant variables are functions of others.In doseresponse studies, Y represents the time to occurrence of a toxic response and x represents the dose metameter, the hazard can be expressed in the form (2) with x Prentice et al. [12] gave specific applications of the Cox model to the analysis of dose-response experiments.The detailed account of dose-response models is available in an expository paper by Kalbfleisch et al. [13].
Bayesian and hierarchical Bayesian estimation of the parameters β 0 , β 1 , . . ., β p , the hazard rate, and the survival function are presented here under the assumptions of the squared error loss function (SELF) and suitable joint prior density of (β 0 , β 1 , . . ., β p ).A numerical illustration based on the model ( 2) is presented for survival data set on advanced lung cancer patients.

Model Parameters
Under the model assumptions ( 1) and ( 2), the hazard rate (HR) and the survival function (SF) are respectively given by The SF (5) gives the probability of survival of an individual with a given vector x  of concomitant variables, up to time t measured from the chosen origin, which may be the start of the clinical study or the point at diagnosis.

Data Set
It is assumed that 'n' individuals enter the clinical study at different points of time and the clinical study lasts a predetermined follow-up period t = T 0 .Let 'd' be the number of individuals responding prior to the follow-up period T 0 , then the rest of individuals, say s = (n -d) consist of those who are lost to follow-up at different time points during the study and those who did not respond till the end of the clinical study.This type of censoring is also known as "progressive censoring" in the literature.It is also assumed that measurements on (p + 1) concomitant variables for all the patients are also available.For this situation the sample data will consist of the observation vectors (t j , x j0 , x j1 , . . ., x jp ), j = 1, 2, . . ., d and (t ⁄ k , x ⁄ k0 , x ⁄ k1 , . . .x ⁄ kp ), k = 1, 2, . . ., s, where t j denotes the time-to-response of the j th individual measured from his entry point, t ⁄ k denotes the censored response-time of the k th individual and x jr and x ⁄ kr , r = 0, 1, 2, . . ., p denote r th concomitant variable on the said j th and k th individuals respectively.

Likelihood Function
For the hazard model ( 2) and the Type III censored sam-ple data set described earlier LF works out as where The product term in ( 6) can be written as a sum as where ∑* is the sum over all possible combinations of m 0 , m 1 , . . ., m p , such that and for given (m 0 , m 1 , . . ., m p ), m S  , d has been defined in the appendix.Hence, the LF can be written as Throughout this paper, g and g* will be used as the generic notations for the prior and the posterior densities respectively and the loss structure will be characterized by the usual squared error loss function (SELF).We shall also use the generic notation K for the normalization constant.

BAYESIAN ESTIMATION OF THE MODEL PARAMETERS
Here it is assumed that prior densities of β 0 , β 1 , . . ., β p mentioned earlier are a priori independent and that β r , r = 0, 1, 2, . . ., p, follows the gamma prior density with known scale and shape hyperparameters b r and a r respectively.For this situation, the joint prior density of (β 0 , β 1 , . . ., β p ) is given by The Bayesian results for the non-informative prior quasi-density (NPQD) specified by are also obtained.The role of NPQD in Bayesian analysis is elucidated in a basic paper of Bhattacharya [14].The raison d ⁄ ệtre for priors mentioned above and details of their use are available in Raiffa and Schlaifer [15].
On the basis of sample data set described earlier, the Bayesian posterior density of ( 0 ,  1 , . . .,  p <∞) is obtained by combining the LF (10) and joint prior density (11) with the help of the Bayes theorem.This works out to be where From ( 13), the marginal posterior density of  r (r = 0, 1, . ., p) is given by Under the assumption of the SELF, the Bayes estimator of  r (r = 0, 1, 2, . .., p) and its posterior variance respectively are obtainable from the following expressions The Bayes estimator of the SF is given by the expression: Evaluating the above multiple integral we obtain Similarly, the posterior variance of Ŝ (t; x) is obtained as The Bayes estimator of the HR of the patient having concomitant variable vector x = x 0 , x 1 , . . ., x p ) is obtained as p r r r 0 where r β  (r = 0, 1, 2, . . ., p) is to be substituted from (16).The Bayesian results for the NPQD ( 12) can be obtained from those corresponding to the prior density (11) given above, by replacing b r = 0, and a r = 1, r = 0, 1, 2, . . .p.

HIERARCHICAL BAYESIAN ESTIMATION OF THE MODEL PARAMETERS
In hierachical Bayes approach [16][17][18][19][20] a second stage prior is assumed for the unknown hyperparameter of the prior distribution assumed at the first stage.Here the hirachical Bayes (HB) estimators of the parameters:  r (r = 0, 1, 2, . . ., p) , the HR and the SF have been derived under the assumptions of the SELF and the gamma distributions as prior and hyperprior densities.At the first stage  r (r = 0, 1, 2, . . ., p) is assumed to follow the gamma prior density.

     
For the unknown scale hyperparameter b r in (22) the following GPD as hyperprior is assumed at the second stage From ( 22) and ( 23) the prior pdf of  r is obtained as Assuming a priori independence, the joint prior density of ( 0 ,  1 , . . .,  p ) in this case comes out to be On the basis of sample data set described earlier obtained by combining the LF (10) joint prior density (25) with the help of the Bayes theorem.This works out to be where K can be computed by using the result (A.2) of the appendix as where From ( 26) the marginal posterior density of  r (r = 0, 1, . ., p) is obtained as where Using the BPD (30) and the result (A.3) of the appendix, the HB estimator of  r (r = 0, 1, 2, . . ., p) is given by Similarly, the posterior variance of r β  (r = 0, 1, 2, . . ., p) is computed as The HB estimator of S and its posterior variance respectively are evaluated as and The HB estimator of the HR for the patient having concomitant variable vector x = (x 0 , x 1 , . . ., x p ) is given by p r r r 0 where r β  (r=0, 1, 2, . . ., p) is to be substituted from (34).

NUMERICAL ILLUSTRATION ON LUNG CANCER PATIENTS
To illustrate the use of the model characterized by the Eq.4, the survival data set on 137 advanced lung cancer patients previously studied by Prentice [21] is used.Patients were randomized according to one of two chemotherapeutic agents: standard and test.To study the possible differential effects of therapy on tumor cell type, tumors were classified into four broad groups termed as squamous, small, adeno and large.The author used four covariates: performance status, time from diagnosis to study, age, and previous therapy to be denoted by x 1 , x 2 , x 3 and x 4 respectively.Assuming the intercept to be zero and taking noninformative prior quasi-density (NPQD) at (12)     tumor types for the standard and test therapies are shown in Table 1.The estimates of the survival function S = S (t; x) for given vector x   From the comparison of estimates for the two therapies, for given arbitrary concomitant vector x  for squamous and adeno tumor cell types, the test therapy prolong the survival of patients.The test therapy comes out to be the most effective for adeno tumor cell type for this particular case.

Figure 1 .
Figure 1.Estimates of the survival functions of standard and test therapies for 'squamous' tumor cell.

Figure 2 .
Figure 2. Estimates of the survival functions of standard and test therapies for 'small' tumor cell type.

Figure 3 .
Figure 3. Estimates of the survival functions of standard and test therapies for 'adeno' tumor cell type.

Figure 4 .
Figure 4. Estimates of the survival functions of standard and test therapies for 'large' tumor cell type.comitant variables are presented inTable 2. Figures 1 to 4 provide the plots of estimates of the survival function of standard and test therapies for different tumor cell types.

Table 1 .
for ( 1 ,  2 ,  3 ,  4 ) the Bayes estimates of the model parameters are obtained for different tumor types and the results for the standard and the test therapies are compared.The Bayes estimates of  1 ,  2 ,  3 , and  4 for different Estimates of model parameters for different tumor types.

Table 2 .
Estimates of the survival function for different tumor types.

Table 2 .
Figures 1 to 4 provide the plots of estimates of the survival function of standard and test therapies for different tumor cell types.