Paper Menu >>
Journal Menu >>
J. Biomedical Science and Engineering, 2010, 3, 711-718 JBiSE doi:10.4236/jbise.2010.37095 Published Online July 2010 (http://www.SciRP.org/journal/jbise/). Published Online July 2010 in SciRes. http://www.scirp.org/journal/jbise Bayesian and hierarchical Bayesian analysis of response - time data with concomitant variables Dinesh Kumar Department of Community Medicine, Government Medical College, Chandigarh, India. Email: dinesh_walia@rediffmai.com Received 9 January 2010; 19 January 2010; 30 January 2010. ABSTRACT This paper considers the Bayes and hierarchical Bayes approaches for analyzing clinical data on re- sponse times with available values for one or more concomitant variables. Response times are assumed to follow simple exponential distributions, with a dif- ferent 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 dose- response modeling have also been explored. Bayesian estimators derived in this paper are applied to lung cancer data set with concomitant variables. Keywords: Bayes Estimator; Bayesian Posterior Density; Gamma Prior Density (GPD); Hierarchical Bayes Esti- mator; Hyperprior; Noninformative Prior Quasi-Density (NPQD); Progressive Censoring; Squared Error Loss Function (SELF); Whittaker Function W s1, s2 (.). 1. INTRODUCTION In biomedical studies, a considerable interest is laid upon developing statistical techniques for analyzing sur- vival data which utilize information available on con- comitant variables. In classical analysis of survival data, several models [1-7] are used for such situations. The usual proportional hazards (PH) regression model pro- posed by Cox [8] has been extensively discussed in the literature. Byar et al. [9] and Greenberg et al. [10] pre- sented 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 ad- justing the effect of a single concomitant variable in the Bayesian framework. The present paper presents the Bayesian and hierarchical Bayesian analysis of re- sponse-time data in more general situations of more than one concomitant variables available for their ef- fects to be adjusted. The exponential survival model -λy fy | λλe0; λ0, (1) representing the death density function (DDF) corre- sponding to the survival time Y is assumed. We also as- sume that the hazard λ for a patient under clinical inves- tigation is linearly related to measurements on ‘p’ con- comitant variables x1, x2, . . . , xp as follows pp 0rr rr01p r=1 r=1 λλt; xββxβx0β,β,...,β (2) where β0, β1, . . . . , βp are (p + 1) unknown parameters and x0 = 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 inter- cept. Of course, it is necessary that the right hand side of (2) be positive. The above hazard model can also be written as λt;xx β (3) where x = (x0 , x1 , . . . . , xp) is a 1 × (p + 1) vector of concomitant variables measured on the individual under clinical investigation and β = (β 0, β1, . . . , βp) ⁄ is a (p + 1) × 1 vector of unknown parameters. A natural extension of the model (2) is a dose-response model with hazard as a polynomial function of con- comitant variables covering the situations wherein some concomitant variables are functions of others. In dose- response 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 = (x0, x1, . . . , xq) and β = (β 0, β1, β2. . . , βq) ⁄ , where q is the number of stages in the dose-response phenomena. D. Kumar / J. Biomedical Science and Engineering 3 (2010) 711-718 Copyright © 2010 SciRes. JBiSE 712 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 avail- able 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 sur- vival 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 illus- tration based on the model (2) is presented for survival data set on advanced lung cancer patients. 2. TOOLS AND TECHNIQUES 2.1. Model Parameters Under the model assumptions (1) and (2), the hazard rate (HR) and the survival function (SF) are respectively given by p rr r=0 ht;x βx (4) p rr r=0 SSt;x exp-βxt t0 (5) The SF (5) gives the probability of survival of an in- dividual with a given vector x of concomitant vari- ables, up to time t measured from the chosen origin, which may be the start of the clinical study or the point at diagnosis. 2.2. 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 = T0. Let ‘d’ be the number of individuals responding prior to the follow-up period T0, 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 re- spond till the end of the clinical study. This type of cen- soring 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 avail- able. For this situation the sample data will consist of the observation vectors (tj, xj0, xj1, . . . , xjp), j = 1, 2, . . . , d and (t ⁄ k, x⁄ k0, x⁄ k1, . . . x⁄ kp), k = 1, 2, . . . , s, where tj de- notes the time-to-response of the jth individual measured from his entry point, t⁄ k denotes the censored re- sponse-time of the kth individual and xjr and x⁄ kr, r = 0, 1, 2, . . . , p denote rth concomitant variable on the said jth and kth individuals respectively. 2.3. Likelihood Function For the hazard model (2) and the Type III censored sam- ple data set described earlier LF works out as p rr r=0 -βQ d 0j0 1j1pjp j=1 01 p βΠβx+ βx... βxe 0β,β,...,β (6) where ds/k /kr rjjr j=1k 1 Qt xtx,r0,1,2,...,p (7) The product term in (6) can be written as a sum as p 01 dm mmp 0j01j1pjp0 1 m,d j=1 Πβx+ βx...βx*Sββ...β (8) where ∑* is the sum over all possible combinations of m0, m1, . . . , mp, such that p rr r=0 mdm0,1,2,...., d (9) and for given (m0, m1, . . . , mp), m S , d has been defined in the appendix. Hence, the LF can be written as rrr pm-βQ r01p m,d r0 β*S βe0β,β,. ..,β (10) 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 normaliza- tion constant. 3. 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 br and ar re- spectively. For this situation, the joint prior density of (β0, β1, . . . , βp) is given by rrr pa1-bβ 01 pr01p r0 gβ,β,...,ββe0β,β,...,β (11) The Bayesian results for the non-informative prior quasi-density (NPQD) specified by 01 p01p gβ,β, ...,β10β,β,...,β (12) are also obtained. The role of NPQD in Bayesian analy- sis is elucidated in a basic paper of Bhattacharya [14]. The raison d ⁄ệtre for priors mentioned above and details D. Kumar / J. Biomedical Science and Engineering 3 (2010) 711-718 Copyright © 2010 SciRes. JBiSE 713 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 ob- tained by combining the LF (10) and joint prior density (11) with the help of the Bayes theorem. This works out to be rr r rr p-βbQ am-1 -1 01 pr m,d r0 01 p g* β,β,...,βK*Sβe 0β,β,...,β (13) where rr prr am m,d r0 rr am K*S bQ (14) From (13), the marginal posterior density of r (r = 0, 1, . . , p) is given by rr r rr uu r p -βbQ uu am-1 -1 ram m,d ru0 uu r g* β am K*Sβe bQ 0β; r=0, 1, 2, ..., p (15) Under the assumption of the SELF, the Bayes estima- tor of r (r = 0, 1, 2, . . ., p) and its posterior variance respectively are obtainable from the following expres- sions uu puu rr -1 ram m,d u0 rr uu am am βK*SbQ bQ (16) uu rrrr -1 r2 m,d rr 2 puu r am u0 uu amam1 VβK*S bQ am β bQ (17) The Bayes estimator of the SF is given by the expres- sion: 00 11pp 00 0 01p0 1p St;x.. exp-βxβx... βxt g* β,β,...,βdβdβ..dβ (18) Evaluating the above multiple integral we obtain rr prr -1 am m,d u0 rrr am St;x K*S bQ+xt (19) Similarly, the posterior variance of Ŝ (t; x) is obtained as rr 2 prr -1 am m,d u0 rrr VS t;x am K*S -St;x bQ+2xt (20) The Bayes estimator of the HR of the patient having concomitant variable vector x = x0, x1, . . . , xp ) is ob- tained as p rr r0 hxβ (21) 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 br = 0, and ar = 1, r = 0, 1, 2, . . . p. 4. HIERARCHICAL BAYESIAN ESTIMATION OF THE MODEL PARAMETERS In hierachical Bayes approach [16-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 dis- tributions as prior and hyperprior densities. At the first stage r (r = 0, 1, 2, . . . , p) is assumed to follow the gamma prior density. r rrr a a1 bβ r rrr r r rrrr b gβgβ|b βe a 0β; b ,a0; a known (22) For the unknown scale hyperparameter br in (22) the following GPD as hyperprior is assumed at the second stage r rrr v v1 ρ b r rrrr r r rrrrr ρ gbgb|p, vbe v 0b; ρ,v0; ρ,v known (23) From (22) and (23) the prior pdf of r is obtained as rrrrrrrr 00 gβgβ, bdbgβ|bgbdb (24) Assuming a priori independence, the joint prior den- sity of (0, 1, . . . , p) in this case comes out to be r r a1 p r 01 pρ r0 rr β gβ, β, ..., βρ+β (25) On the basis of sample data set described earlier ob- tained by combining the LF (10) joint prior density (25) D. Kumar / J. Biomedical Science and Engineering 3 (2010) 711-718 Copyright © 2010 SciRes. JBiSE 714 with the help of the Bayes theorem. This works out to be rr rr r ma-1-βQ p -1 r 01 pv m,dr0 rr βe g* β, β, ..., βK*S ρ+β (26) where K can be computed by using the result (A.2) of the appendix as rr 1 rr 21 1 ~ ρQ W δ-1 2 ρQ2 prr r m,d δ+1 2 r=0 r -δ δ Γm+a eρ K= *S, 22 Q (27) where 1rrr δmav (28) and 2rrr δ1v m a (29) From (26) the marginal posterior density of r (r = 0, 1, . . , p) is obtained as rr r r ~ βQ mrar 1 -1 rv m,d rr βe g* βK*SA(m) ρβ (30) where uu 1 uu 21 1 ρQ δ*12 ρQ2 dδ* uu u δ12 ru=0 u A( m ) -δ* W mae ρ, 22 Q (31) where 1uuu δ*mav (32) and 2uuu δ*1vma (33) Using the BPD (30) and the result (A.3) of the appen- dix, the HB estimator of r(r = 0, 1, 2, . . . , p) is given by rr 2-1 1 rm,drr 1/ 2ρQ δ1 r r βK*SB(m)m+a W-δ1 ρ, Q2 2 (34) where 1 rr 1 δ-12 ρQ2 rr r δ+1 2 r Γm+aeρ B(m) A(m) Q (35) Similarly, the posterior variance of r β (r = 0, 1, 2, . . . , p) is computed as ~ rr 2 1 rrrrr m,d 1/2 ρQ δ21 r r VβK*SB(m)m+am+a1 W-δ2 ρ, Q2 2 (36) The HB estimator of S and its posterior variance re- spectively are evaluated as rr r rr r1 21 1 ~ ρQxt ρQxtδ-1 2 pδ rr r m,d δ+1 2 r=0 rr S(t; x) -δ W Γm+a eρ *S , 22 Qxt (37) and ~ rr r rr r1 21 1 m,d ρQ2xt ρQ2xtδ-12 pδ rr r δ+1 2 r=0 rr 2 V(St; x)*S -δ W Γm+aeρ, 22 Q2xt S(t; x) (38) The HB estimator of the HR for the patient having concomitant variable vector x = (x0, x1, . . . , xp) is given by p rr r0 h=x β (39) where r β (r=0, 1, 2, . . . , p) is to be substituted from (34). 5. 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. Pa- tients were randomized according to one of two chemo- therapeutic agents: standard and test. To study the possi- ble 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 x1, x2, x3 and x4 respectively. Assuming the intercept to be zero and taking noninformative prior quasi-density (NPQD) at (12) 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 D. Kumar / J. Biomedical Science and Engineering 3 (2010) 711-718 Copyright © 2010 SciRes. JBiSE 715 Table 1. Estimates of model parameters for different tumor types. Estimates of i Type of Tumor Type of therapy i Squamous Small Adeno Large 1 2.331 × 10-5 1.841 × 10-5 3.707 × 10-5 1.027 × 10-5 2 4.274 × 10-5 6.649 × 10-5 6.690 × 10-5 1.010 × 10-4 Standard 3 2.206 × 10-5 1.484 × 10-5 6.316 × 10-5 1.520 × 10-5 4 2.091 × 10-4 2.364 × 10-4 8.790 × 10-3 1.802 × 10-4 1 5.753 × 10-6 3.352 × 10-5 3.539 × 10-5 1.739 × 10-5 2 5.041 × 10-5 1.445 × 10-4 2.934 × 10-4 1.561 × 10-4 Test 3 7.354 × 10-6 6.449 × 10-5 3.584 × 10-5 2.116 × 10-5 4 6.563 × 10-5 8.231 × 10-4 7.890 × 10-4 2.201 × 10-4 Table 2. Estimates of the survival function for different tumor types. Estimates of the SF Type of Tumor Therapy Survival time in days Squamous Small Adeno Large Standard 10 0.9488 0.9515 0.4405 0.9583 20 0.9007 0.9057 0.2426 0.9186 30 0.8554 0.8627 0.1506 0.8810 40 0.8128 0.8221 0.1008 0.8451 50 0.7726 0.7837 0.0709 0.8109 60 0.7347 0.7475 0.0518 0.7784 70 0.6990 0.7133 0.0389 0.7475 80 0.6653 0.6810 0.0298 0.7179 90 0.6335 0.6504 0.0234 0.6898 100 0.6034 0.6215 0.0185 0.6630 Test 10 0.9810 0.8580 0.8637 0.9422 20 0.9625 0.7402 0.7497 0.8881 30 0.9443 0.6417 0.6536 0.8376 40 0.9265 0.5589 0.5722 0.7904 50 0.9091 0.4887 0.5028 0.7462 60 0.8921 0.4290 0.4434 0.7048 70 0.8755 0.3779 0.3923 0.6660 80 0.8592 0.3339 0.3482 0.6297 90 0.8432 0.2960 0.3099 0.5956 100 0.8276 0.2631 0.2765 0.5636 D. Kumar / J. Biomedical Science and Engineering 3 (2010) 711-718 Copyright © 2010 SciRes. JBiSE 716 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10 20 304050 60 70 80 90100 Survival time (Days) Estimates of the SF Test therapy Standard therapy Figure 1. Estimates of the survival functions of standard and test therapies for ‘squamous’ tumor cell. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10 2030 405060 70 8090100 Survival time (Days) Estimates of the SF Test therapy Standard therapy Figure 2. Estimates of the survival functions of standard and test therapies for ‘small’ tumor cell type. 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 = (60, 9, 63, 10) (say) of con- 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1020 30 4050607080 90100 Survival time (Days) Estimates of the SF Test therapy Standard therapy Figure 3. Estimates of the survival functions of standard and test therapies for ‘adeno’ tumor cell type. 0 0.2 0.4 0.6 0.8 1 1.2 10 203040 50 607080 90100 Survival time (Days) Estimates of the SF Test therapy Standard therapy Figure 4. Estimates of the survival functions of standard and test therapies for ‘large’ tumor cell type. comitant variables are presented in 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. D. Kumar / J. Biomedical Science and Engineering 3 (2010) 711-718 Copyright © 2010 SciRes. JBiSE 717 From the comparison of estimates for the two thera- pies, 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. 6. ACKNOWLEDGEMENTS The present work is respectful homage to my respected teacher Late Professor Sameer K Bhattacharya, Head, Department of Mathematics and Statistics, University of Allahabad India. The paper was written when author was working asAssistant Professor, Department of Social and Preventive Medicine, MLN Medical College, Allahabad. The help rendered by Dr. Nand Kishore Singh, Department of Mathematics and Statistics, University of Allahabad, at each stage of the present work is greatly acknowledged. The work in the present form was also not possible without critical evaluation and computer programming by Dr. AK Lal, Indian Institute of Technology, Kanpur India. The author is also thankful to Mr. Parminder Kumar, Data Entry Operator, GMCH-32, Chandigarh for his help in typing the manuscript. REFERENCES [1] Cox, D.R. (1964) Some applications of exponential or- dered scores. JRSS, Series B, 26, 103-110. [2] Feigl, P. and Zelen, M. (1965) Estimation of exponential survival probabilities with concomitant information. Biometrics, 21(4), 826-838. [3] Zippin, C. and Armitage, P. (1966) Use of concomitant variables and incomplete survival information in the es- timation of an exponential survival parameter. Biometrics, 22, 665-672. [4] Glasser, M. (1967) Exponential survival with covariance. Journal of American Statistical Association, 62, 561-568. [5] Cox, D.R. and Snell, E.J. (1968) A general definition of residuals, (with Discussion). Journal of the Royal Statis- tical Society B, 30, 248-275. [6] Crowley, J. and Hu, M. (1977) Covariance analysis of heart transplants survival data. Journal of American Sta- tistical Association, 72(357), 27-36. [7] Mantel, N. and Myers, M. (1971) Problems of conver- gence of maximum likelihood iterative procedures in multiparameter situation. Journal of American Statistical Association, 66(335), 484-491. [8] Cox, D.R. (1972) Regression models and life tables, (with discussion), Journal of the Royal Statistical Society B, 34(2), 187-220. [9] Byar, D.P., Huse, R., Bailar, J.C. and The Veterans Ad- ministration Cooperative Urological Research Group (1974) An exponential model relating censored survival data and concomitant information for prostatic cancer patients. Journal of National Cancer Institute, 52(2), 321-326. [10] Greenberg, R.A. and Bayard, S. (1974) Selecting con- comitant variables using a likelihood ratio step-down procedure and a method of testing goodness of fit in an exponential survival model. Biometrics, 30(4), 601-608. [11] Bhattacharya, S.K., Singh, N.K. and Kumar, D. (1995) Bayesian survival analysis of clinical data using a co- variate. Revista Brasileira de Probabilidade e Estadistica, 9, 141-156. [12] Prentice, R.L., Peterson, A.V. and Marek, P. (1982) Dose mortality relationships in RFM mice following 137 ga- mma ray irradiation. Radiation Research, 90(1), 57-76. [13] Kalbfleisch, J.D., Krewski, D.R. and Van Ryzin, J. (1983) Dose-response models for time to response toxicity data. The Canadian Journal of Statistics, 11(1), 25-49. [14] Bhattacharya, S.K. (1967) Bayesian approach to life testing and reliability estimation. Journal of American Statistical Association, 62, 48-62. [15] Raiffa, H. and Schlaifer, R. (1961) Applied statistical decision theory. Harward University Press, Boston. [16] Berger, J. (1986) Multivariate estimation: Bayes, empiri- cal bayes and stein approaches. Society for Industrial and Applied Mathematics, Philadelphia. [17] Bhattacharya, S.K. and Tiwari, R.C. (1992) Hierarchical Bayesian reliability analysis using Erlang families of priors and hyperpriours. Microelectronics Reliability, 32(1-2), 241-247. [18] Bhattacharya, S.K., Singh, N.K. and Tiwari, R.C. (1992) Hierarchical Bayesian survival Bayesian survival analy- sis based on generalized exponential model. Metron, 50(3), 161-183. [19] Good, I.J. (1980) Some history of the hierarchical bayes- ian methodology in bayesian statistics. Bernardo, J.M., De Groot, M.H., Lindley, D.V. and Smith, A.F.M., Eds., University Press, Valencia. [20] Stangl, D.K. and Greenhouse, J.B. (1998) Assessing placebo response using Bayesian hierarchical survival models. Lifetime Data Analysis, 4(1), 5-28. [21] Prentice, R.L. (1973) Exponential survival with censor- ing and explanatory variables. Biometrika, 60(2), 279- 288. [22] Gradshteyn, I.S. and Ryzhik, I.M. (1980) Tables of Inte- grals, Series and Products. Academic Press, New York. D. Kumar / J. Biomedical Science and Engineering 3 (2010) 711-718 Copyright © 2010 SciRes. JBiSE 718 APPENDIX The definition of m,d S notation and a mathematical result on a special function used in the present work are presented here. 1) Them,d S Notation Let X denotes a d x (p + 1) matrix of concomitant variables given below 10 111p 20 212p d0 d1dp xx...x xx...x X = xx...x (A.1) and m = (m0, m1, . . . . , mp) . For given combinations of m0, m1, . . . . , mp, satisfying the condition: m0 + m1 + . . . + mp = d (A.2) m,d S is given by the sum of all products of mr ele- ments from rth column (r = 0, 1, 2, . . . , p) of the matrix X, such that no two elements in the product term lie in same row of the matrix. 2) Integral for the Whittaker Function The following variant of the integral representation ([22] p.319, Sec. 3. 383, Formula 4) has been used in the present work: 1-q-A q-A , 22 A-1 -Bw q-A-1/2A-q-1 /2 Bv/2 W q 0 we dw = A eBvBv v+w (A.3) This results holds good provided that ReA > 0, ReB > 0, where 12 s, s W( . ) is the well known Whittaker func- tion. |