Inference for a General Class of Models for Recurrent Events with application to cancer data

Inference for a General Class of Models for Recurrent Events with application to cancer data Juan R. Gonz´alez Ruiz Center for Genomic Regulation and Catalan Institute of Oncology Academic Dissertation for the Degree of PhD advised by Dr. Edsel A. Pe˜na The University of South Carolina at Columbia and Dr. Pedro Delicado Useros Universitat Polit`ecnica de Catalunya PhD thesis submitted to Universitat Polit`ecnica de Catalunya Barcelona, April 2006 To Manoli, Iker and Jan To my parents Acknowledgments I specially wish to express my deepest gratitude to Professor Edsel Pe˜na for his guidance and support all through the years of working on this PhD thesis. In spite of having a lot of work he spent a lot of his time reading hundreds of e-mails I sent him. He always had kind words about the work I did and encouraged me to continue working hard on my research. I am indebted to him, not only for everything I have previously mentioned, but also for being an excellent mentor, introducing me as a scientist, encouraging me to travel and bringing me into contact with many important researchers. I also want to thank Professor Pedro Delicado for his guidance and for co-supervising this thesis with such great interest. He has supported me with advices and ideas, and he was always welcome to discuss my work and future directions with me. A very special thanks goes out to Dr. Merc`e Peris, head of Cancer Prevention and Control Unit from Catalan Institute of Oncology. She encouraged me to work on recurrent events settings applied in cancer data. She was very kind allowing me to make my work compatible in her unit with this PhD thesis. I also would like to thank the two anonymous reviewers for the careful reading of the manuscript and for providing me lots of invaluable comments, constructive criticisms, and very interesting suggestions. Their helpful and thoughtful comments have improved many aspects of this thesis dissertation. I like to thank Dr. Elizabeth Slate for her helpful feed-back on some parts of this dissertation. I would specially like to acknowledge her comments about the dynamic cancer model. I am also very grateful to her because she sent me some very useful Rfunctions to summarize simulations iii iv Chapter 0 Acknowledgments results showed in Chapter 5. We thank Dr Octavi Servitge and Dr Cristina Muniesa for providing the non-Hodgkin’s lymphoma data set and for their special effort in obtaining the treatment response information. I must also acknowledge Professor Narc´ıs Nabona for his helpful conversations on maximization problems, particularly on the topic about penalized approach. He was very kind to supply me with some Fortran routines to solve optimization problems. Very special thanks I address to Dr. Virginie Rondeau with whom I shared ideas about penalized likelihood estimation for correlated survival data. She very kindly sent me the Marquart algorithm programmed in Fortran. I would also like to express my gratitude to all people at the Cancer Prevention and Control Unit from Catalan Institute of Oncology. I would specially like to mention Merc`e Margalef who always encouraged me to continue on my graduate career. US NSF Grant DMS 02080301 is acknowledge because it allowed me to attend the International Conference on Reliability and Survival Analysis 2003 and met Professor Pe˜na. Grant to attend the seminars on Empirical processes: theory and statistical applications by UNESCO, held in EMS Summer School of Laredo is also acknowledge. Es bastante dif´ıcil poder agradecer con un ´unico p´arrafo la ayuda que Xavier Sol´e me ha prestado a lo largo de estos ´ultimos a˜nos. Quiz´as su ayuda con el ingl´es haya sido una mera an´ecdota con el resto de cosas que he aprendido con ´el. Recuerdo que empez´o ayudandome con el “maravillos” mundo de la creaci´on de dlls hace m´as de cuatro a˜nos. Perdimos muchas horas hasta que conseguimos llamar por primera vez a Fortran desde R. Tambi´en he recibido mucha ayuda a la hora de programar algunas cosas pues sus conocimientos en inform´atica (a parte de otros campos de la ciencia) no parecen tener l´ımite. Pero sin duda, su mejor ayuda a esta tesis, ha sido el estar disponible en cualquier momento. He podido hablar con ´el cuando necesitaba a alguien y siempre me ha dado excelentes consejos sobre todo cuando mi ´animo deca´ıa. Encara que ja no estiguem treballant junts, ja saps que mai no deixarem de veure’ns. Persones com tu s´on les que realment v fan que la paraula AMIC tingui sentit i no sigui m´es que una definici´o al diccionari. Espero que la vida et torni tot all`o que tu dones i que no canvi¨ ıs mai. No te n’oblidis que tenim una cita anual al Xin´es!. Tambi´en quisiera agradecer especialmente la ayuda que he recibido de Elisabet (Looking). Ya he perdido la cuenta de los p´arrafos en ingl´es que me ha corregido. Tambi´en me gustar´ıa agradecerle su amistad i as´ı como las numerosas charlas que hemos tenido. Entre ellas, ha tenido que soportar algunos de mis rollos estad´ısticos. A Klaus Langohr quiero agradecer su enorme paciencia que ha tenido al ayudarme a escribir esta tesis en LaTeX. Realmente ha sido como un manual on-line, en el buen sentido de la palabra. Siempre ha estado dispuesto a solucionar cualquier duda que me ha surgido. De hecho, todo lo poco que s´e de LaTeX me lo ha ense˜nado ´el. A pesar que s´olo llevo unos pocos meses trabajando en el Centro de Regulaci´on Gen´omica, tambi´en hay personas trabajando conmingo que merecen estar en estos agradecimientos. Llu´ıs Armengol y M´onica Gratac´os quiz´as sean las personas que m´as me han soportado durante los ´ultimos meses en la realizaci´on de esta tesis. Este periodo final coincide sin duda con los meses m´as duros pues es cuando hay que preparar el manuscrito final y su defensa. Ambos me han ayudado a superar esta parte final de la tesis anim´andome y sobretodo d´andome excelentes consejos. Finalmente, y no por ser lo ´ultimo es lo que menos ganas tengo de agradecer, quisiera mostrar mi mas sincero agradecimiento a mis padres. Desde que mi memoria puede recordar s´olo veo el esfuerzo que realizaron para que todos sus hijos pudieran tener una educaci´on que ellos no pudieron recibir. Mis ojos se llenan de l´agrimas al recordar a mi padre queriendo que estudiara ya fuera verano o invierno, hubieran o no vacaciones. Este momento es el que culmina su prop´osito. Esta tesis doctoral est´a dedicada completamente a mi mujer, Manoli, que me ha apoyado durante estos ´ultimos a˜nos en todo momento y que siempre me ha animado a realizar estos estudios. Su amor, cari˜no, y ternura han hecho que este periodo de tiempo sea como un paseo por un precioso bosque lleno de bonitos lugares para tomar aliento. Nuestros hijos, Iker y Jan, han sido fruto del lugar m´as bello, sin duda. Abstract Survival analysis arises when we are interested in studying statistical properties of a variable which describes the time to a single event. This type of analysis occurs commonly in two areas: biomedicine and engineering. In biomedicine research it is known as survival analysis and refers often to the time from the beginning of the treatment to the occurrence of a particular condition or death. In some situations, we may observe that the event of interest occurs repeatedly in the same individual, such as when a patient diagnosed with cancer tends to relapse over time or when a person is repeatedly readmitted in a hospital. In this case we speak about survival analysis with recurrent events. Recurrent nature of events makes necessary to use other techniques from those used when we analyze survival times from one single event. In this dissertation we deal with this type of analysis mainly motivated by two studies on cancer research that were created specially for this research. One of them belongs to a study on hospital readmissions in patients diagnosed with colorectal cancer, while the other one deals with patients diagnosed with non-Hodgkin’s lymphoma. This last study is mainly relevant since we include information about the effect of treatment after relapses and some authors have stated the needed of developing an specific model for relapsing patients in cancer settings (Montoto et al., 2002). These two data sets together with two other existing examples and that have been extensively analyzed, have been used to illustrate the statistical methodology proposed in this work. In this dissertation, we address two different types of statistical analysis similar to those one may carry out when we deal with survival analysis for a single event. This two types of analysis are known in the biomedical literature as univariate or multivariate analysis. Univariate analysis studies how one variable may modify the probability of observing a new recurrence. Our contribution to this problem is to propose a method to construct confidence intervals for vii Resum La necessitat de l’an`alisi de superviv`encia apareix quan ens cal estudiar propietats estad´ıstiques d’una variable que descriu el temps fins que succeeix un esdeveniment ´unic. Aquest tipus d’an`alisi se sol plantejar normalment en dues `arees: la biomedicina i l’enginyeria. En investigaci´o biom`edica, es coneix com a an`alisi de superviv`encia, i usualment fa esment, al temps des de l’inici del tractament fins a l’esdeveniment d’una condici´o particular o la mort. En algunes ocasions, podem observar que l’esdeveniment d’inter`es se succeeix de manera recurrent en un mateix individu, com pot ser el cas d’un pacient diagnosticat de c`ancer que recau al llarg del temps o quan una persona ´es reingressada repetides vegades en un hospital. En aquest cas, parlem d’an`alisi de superviv`encia amb esdeveniments recurrents. La naturalesa recurrent dels esdeveniments fa necess`aria la utilitzaci´o de t`ecniques estad´ıstiques distintes a les que emprem quan analitzem temps de superviv`encia en el cas d’un esdeveniment ´unic. En aquesta tesi, tractem aquest tipus d’an`alisi motivats per dos estudis de recerca en c`ancer que varen ser dissenyats per dur a terme aquest treball. Un d’ells gira en torn de readmissions hospital`aries en pacients diagnosticats amb c`ancer colorectal, per b´e que l’altre fa refer`encia a pacients diagnosticats amb limfomes no Hogkinians. Aquest darrer estudi ´es especialment rellevant ja que v`arem incloure informaci´o sobre l’efecte del tractament despr´es de les recaigudes, i alguns autors han mostrat la necessitat de desenvolupar un model espec´ıfic per pacients que presenten aquest tipus de malalties (Montoto et al., 2002). Aquest dos conjunts de dades, juntament amb d’altres existents a la literatura biom`edica i que han estat `ampliament analitzats, han estat emprats per il˚ulustrar la metodologia estad´ıstica que es proposa en aquest treball. En aquesta tesi, hem tractat dos tipus diferents d’an`alisi estad´ıstic similars als que es duen a terme quan es tracten dades de superviv`encia simples. Aquests dos tipus d’an`alisi s´on coneguts xv xvi Chapter 0 Resum en la literatura biom`edica com an`alisi univariant i multivariant. L’an`alisi univariant estudia com una variable pot modificar la probabilitat d’observar una nova ocurr`encia. La nostra contribuci´o a aquest problema ha estat proposar un m`etode per construir intervals de confian¸ca per a la mediana de superviv`encia en el cas d’esdeveniments recurrents. ´ Es per aix`o que hem utilitzat dues aproximacions. Una d’elles es basa en les vari`ances asimpt`otiques derivades de dos estimadors de la funci´o de superviv`encia (Pe˜na et al., 2001 y Wang y Chang, 1999), per b´e que l’altre utilitza t`ecniques de remostratge. Aquesta darrera aproximaci´o es ´util ja que un dels estimadors proposats per Pe˜na et al. (2001) encara no t´e una forma tancada per la seva vari`ancia. L’originalitat d’aquest treball rau en el disseny del remostratge quan es tracten dades amb events recurrents, que apareixen amb un esquema conegut com ”sum-quota accrual”i la informativitat del mecanisme de censura per la dreta que presenten aquest tipus de dades. Demostrem la converg`encia d`ebil i constru¨ ım els intervals de confian¸ca assimpt`otics utilitzant aquest resultat. D’altra banda, l’an`alisi multivariant fa esment al problema de com incorporar m´es d’una covariable en l’an`alisi. En problemes amb esdeveniments recurrents tamb´e necessitem tenir en compte que, a m´es de les covariables, l’heterogene¨ ıtat, el nombre d’ocurr`encies, o especialment, l’efecte de les intervencions despr´es de les reocurr`encies poden modificar la probabilitat d’observar un nou esdeveniment en un pacient. Per abordar aquest problema, hem basat el nostre treball en un model proposat per Pe˜na i Hollander (2004). La nostra contribuci´o a aquest punt ha estat l’adaptaci´o de les recaigudes en c`ancer utilitzant aquest nou model per a esdeveniments recurrents en qu`e l’efecte de les intervencions es presenta mitjan¸cant un proc´es anomenat ”edat efectiva”, que actua sobre la funci´o de risc basal. Hem anomenat a aquest model model din`amic de c`ancer (“dynamic cancer model”). Tamb´e tractem el problema de l’estimaci´o de par`ametres de la classe general de models per esdeveniments recurrents proposat per Pe˜na i Hollander (2004), on el model din`amic de c`ancer es pot veure com un cas especial d’aquest model general. Hem desenvolupat dues aproximacions: la primera es basa en infer`encia semiparam`etrica, on la funci´o de risc basal s’especifica de forma no param`etrica; la segona ´es una aproximaci´o basada en versemblan¸ca penalitzada. En el cas d’infer`encia semiparam`etrica adoptem dues estrat`egies diferents depenent de si ajustem un model amb fragilitat (”frailty”) o de si no ho fem. L’algoritme EM ´es utilitzat quan la fragilitat ´es inclosa en el model. Pel que fa a l’aproximaci´o mitjan¸cant versemblan¸ca penalitzada, tamb´e adoptem dues estrat`egies diferents. Una d’elles va ser proposada per Therneau et al. (2003) en el context xvii d’un model de fragilitat compartida (”shared frailty model”). La seva idea es basa en penalitzar la versemblan¸ca parcial en qu`e la penalitzaci´o recau en els coeficients de regressi´o. La segona aproximaci´o basada en penalitzaci´o, tamb´e aplicada en aquest mateix model, va ser proposada per Rondeau et al. (2003). El seu m`etode d’estimaci´o es basa en penalitzar la versemblan¸ca completa i proporciona una estimaci´o no param`etrica de la funci´o de risc basal utilitzant un estimador continu. La soluci´o s’aproxima mitjan¸cant splines. L’avantatge m´es important d’aquest m`etode ´es que podem obtenir f`acilment una estimaci´o suau de la funci´o de risc, aix´ı com una estimaci´o de la vari`ancia de la vari`ancia de la fragilitat, ja que amb les altres aproximacions aix`o no ´es possible. A m´es, aquest darrer m`etode t´e un cost computacional bastant m´es baix que els altres. Les simulacions dutes a terme sota diferents escenaris i mides mostrals, han mostrat bones propietats dels estimadors proposats. A m´es, els resultats obtinguts amb dades reals indiquen que la flexibilitat d’aquest model es una garantia per analitzar dades de pacients que recauen al llarg del temps i que s´on intervinguts despr´es de les recaigudes tumorals. L’aspecte computacional ´es una altra de les contribucions importants d’aquesta tesi al camp dels esdeveniments recurrents. Hem desenvolupat tres paquets de Ranomenats survrec,gcmrec yfrailtypack que s´on accessibles a trav´es de CRAN, http://www.r-project.org/. Aquest paquets permeten a l’usuari calcular la mediana de superviv`encia i els seus intervals de confian¸ca, estimar els par`ametres del model de Pe˜na i Hollander (en particular el model din`amic de c`ancer) mitjan¸cant l’algoritme EM y la versemblan¸ca penalitzada, respectivament. Contents Acknowledgments iii Abstract vii Resumen xi Resum xv 1 Introduction: State-of-the-art 1 1.1 Goals ........................................... 1 1.2 Survival Analysis with Recurrent Events . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Sum-quota accrual scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Doubly-indexed processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Estimation of the survival function . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.1 Effective age process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.2 Reliability models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.3 Pe˜na-Strawderman-Hollander estimator . . . . . . . . . . . . . . . . . . . . 13 1.3.4 Wang-Chang estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.4 Within-subject correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.4.1 Variance-corrected models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4.2 Frailty models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.4.3 Cox extension models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 1.5 General class of models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 xix xx CONTENTS 1.5.1 Pe˜na and Hollander model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.6 ThesisOverview ..................................... 26 2 Studies with Recurrent Event Data 29 2.1 Hospital Readmission Times in Colorectal Cancer . . . . . . . . . . . . . . . . . . . 29 2.1.1 Variables of the data set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.1.2 Descriptive analysis of the data set . . . . . . . . . . . . . . . . . . . . . . . 30 2.2 Non-Hodgkin’s Lymphoma Cancer Relapses . . . . . . . . . . . . . . . . . . . . . . 34 2.2.1 Summary of the data set . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.3 Bladder Cancer Relapses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.4 MMCdataset ...................................... 36 3 Confidence Intervals for Median Survival 39 3.1 Estimation of median survival time and other quantiles . . . . . . . . . . . . . . . 40 3.2 Bootstrapping ζp..................................... 44 3.3 SimulationStudy..................................... 46 3.3.1 Simulation Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4 Examples ......................................... 56 3.4.1 MMCdataset .................................. 56 3.4.2 Colorectal cancer rehospitalizations . . . . . . . . . . . . . . . . . . . . . . . 58 3.4.3 Bladder cancer data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.5 R instructions for survrec package........................... 64 4 Inference for the General Class of Models 71 4.1 Semiparametric inference: EM algorithm . . . . . . . . . . . . . . . . . . . . . . . . 72 4.2 Penalized Likelihood Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.2.1 Penalized partial likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2.2 Penalized full likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.3 Statistical Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 CONTENTS xxi 4.4 Computational Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.5 Hospital Readmission and Bladder Cancer Data Sets Revisited . . . . . . . . . . . 85 4.5.1 Hospital Readmission Study . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.5.2 Bladder Cancer Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.5.3 Miss-specification of effective age . . . . . . . . . . . . . . . . . . . . . . . . 91 4.6 R instructions for gcmrec package ........................... 94 5 Dynamic Cancer Model for Tumor Relapses 97 5.1 The Pe˜na and Hollander Model Revisited . . . . . . . . . . . . . . . . . . . . . . . 99 5.2 Effective Age Process for Cancer Data . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.3 Simulation Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3.1 Simulation Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.3.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.4 The non-Hodgkin’s lymphoma study . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.5 R instructions for gcmrec package ...........................118 6 Concluding Remarks and Future Research 119 6.1 Conclusions........................................119 6.2 FutureResearch .....................................121 Bibliography 122 Appendix 133 A Counting Processes in Survival Analysis 135 A.1 Counting processes approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 A.2 Nonparametric methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 A.3 Cox proportional hazards model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 A.3.1 Stratified Cox model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 xxii CONTENTS B Semiparametric Inference for Pe˜na and Hollander model 145 B.1 Case without Frailties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 B.2 Case with Frailties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 C Published work related to present thesis 151 C.1 Papers...........................................151 C.2 Rpackages.........................................154 C.3 Oral Contributions in Meetings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 D R Functions and Classes 157 D.1 The survrec Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 MMC ...........................................157 Survr ...........................................158 surv.search ........................................159 colon............................................160 mlefrailty.fit........................................161 plot.survfitr........................................164 print.survfitr .......................................165 psh.fit ...........................................166 q.search ..........................................168 summary.survfitr .....................................169 survdiffr..........................................170 survfitr ..........................................172 wc.fit............................................173 D.2 The gcmrec Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 GeneratedData ......................................176 addCenTime .......................................176 gcmrec-internal......................................177 gcmrec...........................................178 CONTENTS xxiii graph.caltimes ......................................183 hydraulic .........................................185 lymphoma.........................................185 plot.gcmrec ........................................186 print.gcmrec........................................187 readmission ........................................188 summary.gcmrec .....................................189 D.3 The frailtypack Package . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 frailtyPenal ........................................191 plot.frailtyPenal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195 print.frailtyPenal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196 readmission ........................................197 summary.frailtyPenal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198 E Additional Rfunctions 201 2Chapter 1 Introduction: State-of-the-art (2002) argued. For this reason, one of the main aims of this dissertation has been to investigate how to incorporate information about interventions in patients who relapse over time. To do so, we used a very flexible model designed for analyzing recurrent event data proposed by Pe˜na and Hollander (2004). In particular, we used the response to the treatment after cancer relapses to model the effect of interventions after reoccurrences. We called this model dynamic cancer model. After that, another important goal was to develop statistical procedures for estimating model parameters involved in this new model such as penalized likelihood methods. Finally, we focused on implementing these methods in a widely used statistical software such as R, to help physicians to carry out their data analysis. In the following sections we outline a survey on statistical methods for the analysis of recurrent event data. In particular, we discuss two important aspects in this type of studies: the heterogeneity and the event dependence. Finally, we illustrate how some models, mainly used in reliability settings, are useful to take into account another important topic which appears in these data: the effect of interventions after re-occurrences. The very flexible model proposed by Pe˜na and Hollander (2004) which is the basis of our work is described at the end of this chapter. 1.2 Survival Analysis with Recurrent Events Survival analysis arises when we are interested in studying statistical properties of the variable T, which describes the time to a single event. This type of analysis occurs commonly in two areas. In medical research it is known as survival analysis and refers often to the time from the beginning of the treatment to the occurrence of a particular condition or death. In engineering it is concerned with reliability and the analysis of failure times. That is, how long a component can be used until it fails. Counting process formulation and martingale theory have become the most used tool in the modern theory of these type of data. In the Appendix A we outline the main points in the use of counting processes as applied in survival analysis. Therneau and Grambsch (2000) can also be looked up for an overview from a very intuitive point of view. However, in many other situations, we observe that the event of interest occurs repeatedly in the same subject such as when a patient diagnosed with cancer tend to relapse over time or when a person is repeatedly readmitted in a hospital. In that case we speak about survival analysis for recurrent events. Repeated events are prevalent in a vast variety of disciplines. These include 1.2 Survival Analysis with Recurrent Events 3 biomedicine, psychiatry, engineering, or sociology among others. Recurrent nature of events makes necessary to use other techniques from those used when we analyze survival times from one single event. In this case, random observation periods for each subject were allowed. Then, we observe the event of interest (always the same) at different times during the follow-up time. Finally, the last event is not observed since it would be observed after the end of study. Thus, censored data appears. This type of scheme which generates the data is known as sum-quota accrual scheme in the literature. This chapter starts by giving a representation of the sum-quota accrual scheme and introducing some notation of this scheme induced by the multiple occurrences and the randomness of follow-up period. The special nature of the data invalidates the direct use of martingales methods. To solve this problem, Section 1.2.2 illustrates another formulation, called doubly-indexed process that was first proposed by Gill (1981) and Sellke (1988) and extended by Pe˜na et al. (2001). This technique will be useful to derive some theory for recurrent events. In Section 1.3 we give a survey of the statistical methodology for recurrent event data, focusing on nonparametric methods. Section 1.4 deals with two important topics which appear in studies with repeated measures: the heterogeneity and the event dependence. Both problems arise when the independent assumption is violated. This section also describes two methods to solve these problems: variance-corrected and frailty models. Finally, Section 1.4.3 describes conditional and marginal models which allow us to incorporate the effect of concomitant covariates, to control the event dependence, and the heterogeneity. 1.2.1 Sum-quota accrual scheme Figure 1.1 illustrates a pictorial representation of our setting, known as sum-quota accrual scheme. We consider a patient diagnosed with cancer (e.g., an observational unit) which is being monitored for the occurrence of a recurrent event over a study period [0, τ] , where τmay represent an administrative time, time of study termination, or some other right-censoring variable (see τin Figure 1.1). The time τcould be a random time governed by an unknown probability distribution function G(t) = Pr(τ≤t). In recurrent event data, times are indexed in two scales, calendar and interoccurrence or gap times. Calendar times are defined by the sequence S0≡0< S1< S2< S3< . . . and correspond to the successive calendar times of event recurrences (see Sin Figure 1.1). Interoccurrence times will be denoted by T1, T2, T3, . . . and correspond to the time between successive event occurrences (see Tin Figure 1.1). Thus, for i= 1,2,3, . . .,Ti=Si−Si−1and 4Chapter 1 Introduction: State-of-the-art Figure 1.1: Graphical representation of the sum-quota accrual scheme for an individual. Si=T1+T2+. . . +Ti.Over the observation period [0, τ], the number of event occurrences is K= max{k∈ {0,1,2, . . .}:Sk≤τ},which is a random variable whose distribution depends on the distributional properties of the inter-occurrence times Tis and the distribution G(w) = P{τi≤w}of τ. As such, Kis informative with regards to the distributional properties of event occurrences (in the example, K= 2). The fact of following-up patient for a fixed time, can lead to some event not being completely observed. Thus, Figure 1.1 illustrates that third event is not observed since it will appear after the end of study. So, T3is not observed completely: we only observe the censored time τ−S2. We notice that the interocurrence times are affected by unobserved variable called frailty (see Zin Figure 1.1). These frailties might make some patients have more recurrences than others depending on their values. In the Appendix A we show that the risk indicator is one of the main components to use in counting processes and martingale theory. We illustrate that, in recurrent event situations, two different time scale are possible. Thus, we have two possible formulations for risk intervals: gap time or total time. Figure 1.2 shows the two types of risk intervals for two hypothetical subjects. For the first patient the interoccurrence times are 5 and 20 months and the third event is not observed after a follow-up of 15 months from the last event (15 is the censored time). These interoccurrence times correspond to calendar times 5, 25 and 40 respectively. For the second patient no events are observed after 30 months. Gap time represents the time from the prior event, e.g. the interoccurrence time. Thus, the gap time formulation (right bottom panel Figure 1.2) for our example indicates that the first subject is at risk of the first event during the interval (0,5], and for the second and third events during (0,20] and (0,15] respectively. On the other hand, total time is the time from a selected point. There exist several examples for this point, among others the date of diagnosis, the time of start of treatment, or the date of birth. In the 1.2 Survival Analysis with Recurrent Events 5 Figure 1.2: Illustrations of the risk interval formulation using an hypothetical data, where ¥ corresponds to and event and •is censoring example (bottom left panel Figure 1.2) , the first subject is at risk for the first, second and third events during the intervals (0,5], (0,25], and (0,40], respectively. Second subject is at risk during (0,30]. 1.2.2 Doubly-indexed processes Counting processes is a powerful tool in survival analysis. However, as Sellke (1988) observed, when dealing with recurrent event data, one should not only consider a calendar time formulation. We also need to consider gap time formulation. These idea originated from Gill (1981) and has been extended in Pe˜na et al. (2001). We begin by defining the following processes, which consider calendar time only. For i= 1, . . . , n and s≥0 let N† i(s) = ∞ X j=1 I(Sij ≤s, Sij ≤τi) (1.1) be processes which count the number of failures for unit iat time swhich have not been censored, and Y† i(s) = I(τi≥s),(1.2) which indicates if unit ihas been censored at time s. Now, we define a filtration G={Gs:s≥0} such that {(N† i(s), Y † i(s)) : s≥0}is G-adapted. Moreover, let A† i(s) = Rs 0Y† i(v)λ(v−SiN† i(v−))dv 6Chapter 1 Introduction: State-of-the-art Figure 1.3: Doubly-indexed processes illustration for an hypothetical case which makes M† i(s) = N† i(s)−A† i(s), s≥0 to be a local square-integrable G-martingale with predictable quadratic covariation process hM† i, N† ii(s) = A† i(s)I{i=i0}(an introduction to these counting processes can be found in Pe˜na et al., 2001). Now, we have the main ingredients to define similar processes to those defined for single event case (see Appendix A) but for recurrent event situations. We introduce appropriate processes that are indexed by calendar time sand gap time t(doubly indexed processes) as follows. This processes are the basic ones considered in Pe˜na et al. (2001) and Sellke (1988) and they provide the connection between the gap time formulation and that based on calendar time. For i= 1,2, . . . , n, let Zi(s, t) = I{s−SiN† i(s−)≤t}be the indicator that for calendar time sat most ttime units have elapsed since the time of the last event. For s, t ≥0 we define Ni(s, t) = Zs 0 Zi(v, t)dN† i(v); Ai(s, t) = Zs 0 Zi(v, t)dA† i(v); (1.3) Mi(s, t) = Zs 0 Zi(v, t)M† i(dv) = Ni(s, t)−Ai(s, t); Yi(s, t) = N† i((s∧τi)−) X j=1 I{Tij ≥t}+I{(s∧τi)−SiN† i((s∧τi)−)≥t}, i = 1,2, . . . , n. We notice that N(s, t) counts the number of observed events occurring over the calendar period [0, s] whose interoccurrence times were at most t. On the other hand, Yi(s, t) counts the number 1.3 Estimation of the survival function 7 of observed events on calendar period [0, s] whose interoccurrence times were at least t. Figure 1.3 shows a hypothetical case followed during 310 months. This patient presents three recurrences at months 110, 185, and 280 from the beginning of study. This fact implies that interoccurrence times are 110, 75, 55, and the censored time correspond to 30 months. Let us assume that we are interested in computing the single processes, N(t) and Y(t) for a selected interoccurrence time t= 100. In this case N(t= 100) = 1 and Y(t= 100) = 2. For the calendar time scale, s= 200, we have N(s= 200) = 2 and Y(t= 200) = 1. Now, let us assume that we would like to know double-indexed processes for both selected interoccurrence and calendar times. Using both time scales we observe that N(s= 200, t = 100) = 1 and Y(s= 200, t = 100) = 1. 1.3 Estimation of the survival function The aim of this section is to give a survey of the statistical methodology for recurrent event data, focusing on nonparametric methods. We will also show how existing reliability models can be applied in biomedical or public health settings. Previously, we need to define a new process called effective age which is outlined in Section 1.3.1. Two estimators for the case of correlated interoccurrence times are described in Sections 1.3.3 and 1.3.4. The methodology described for nonparametric methods in reliability can be found in Hollander and Sethuraman (2002). The estimators for the correlated case are also described and compared in Pe˜na et al. (2001) or in Gonz´alez and Pe˜na (2004). Statistical inference in the presence of recurrent event data has been considered by several authors such as Gill (1981), Vardi (1982a,b), McClean and Devine (1995), Soon and Woodroofe (1996), Wang and Chang (1999) (WC), and Pe˜na, Strawderman, and Hollander (2001) (PSH). A main aspect with this type of data is the sum-quota accrual scheme which leads to an informative stopping rule as well as an informative censoring mechanism (see Section 1.2). Except in PSH (2001), most papers have used restrictive data accrual and censoring schemes for recurrent event data. In PSH (2001) it is assumed that the interoccurrence times represent independent and identically distributed (i.i.d.) observations from an unknown continuous distribution F, and that each subject is observed for a possibly random period of time. As a consequence, the number of event occurrences for a subject or unit is a random variable whose distribution depends on F, hence is informative about it. Moreover, the last observation for each subject is always right- 8Chapter 1 Introduction: State-of-the-art censored, with the censoring variable depending on the length of the observation period and on the previous interoccurrence times for that subject, rendering the censoring mechanism to become informative. 1.3.1 Effective age process Reliability works with models which deal with repairable systems. These models incorporate the idea of effective age, also called virtual age, which describes the effect of repair. Kijima et al. (1988) introduce the notion of the virtual age as an improvement to minimal repair models. Maybe these models are more easily introduced in reliability settings than in biomedical problems due to the ability to monitor the effective age. Perhaps the main reason is that, upon failure, we can decide if the system is replaced by a new identical (perfectly repaired) or if it is restored to its state just before failure (minimally repaired). However, when we are dealing with patients, firstly, we cannot decide the degree of their “repair”: we always want to eliminate all manifestation of the disease (e.g., “perfect” repair); and secondly, we always do not achieve the same response: we can obtain a complete remission of the disease, partial remissions or no response. Although it seems natural to take into account the performed interventions, none of models which are mostly used in biomedical settings, described in Section 1.4.3, incorporate this effect. Before starting to describe how to adopt some existing reliability models in cancer data we focus more on the effect of interventions. For example, people with coronary heart disease are advised to alter their lifestyle by reducing stress level, quitting smoking habit, or doing regular and moderate physical activity. These advises try to modify the probability of presenting a new heart attack. Another example arises from patients with epileptic seizures, where recommendations in order to reduce the number of new seizures include sleeping 8 hours daily, avoiding exposition to flickering lights or reducing alcohol use. Considering cancer settings, in particular some indolent tumors, medical doctors may perform some prophylactic or curative interventions such as chemotherapy, radiotherapy or bone marrow transplantation to improve patient’s disease status. Finally, in reliability the types of interventions are in general simplest since after having a break in the system the piece is repaired or replaced by a new one. There exist several reliability models which incorporate the effect of performed interventions. However, as we have pointed out above, these models have not been applied neither in biomedical nor in public health settings. Next sections give an overview of these models. 1.3 Estimation of the survival function 9 Another justification to use effective age based models is that most of repeated events survival models studied in medicine and public health problems are based on extensions of Cox model (see Section 1.4.3) which assumes that the effect after each intervention is always the same. For example, models which employ total time formulation assumes that all interventions produce a minimal improvement in the patient, e.g. disease continues in a stable manner. Models based on gap time formulation assume that all interventions lead to perfect recovery of the patient, e.g. disease disappears or a complete remission is achieved. However, it is not the case that the effect of treatment intervention will always be the same. For example, patients with a recurrent tumor usually are treated after each relapse and they may obtain a different response after each treatment. Patients sometimes achieve a complete remission, others minimal and others between none and complete remission, e.g. patients suffer a little improvement but the disease still remains (in reliability terms an intervention between perfect and minimal repair). Thus, as we will illustrate in Chapter 5, it can be very important to monitor the effective age process when we are dealing with some biomedical data. Although obtaining information about the effective age in health problems may be complicated, in cancer settings there exist the possibility to get this information using oncological terminology. The adoption of effective age to biomedical data, and in particular to some indolent cancers such as some lymphomas, has been an important part of this thesis. The effective age for the i-th unit is defined as an observable processes {Ei(s):0≤s≤ s∗}, i = 1,2, . . . , n, satisfying the following conditions: (I) Ei(0) = ei0, almost surely (a.s.), where ei0, i = 1,2, . . . , n, are nonnegative real numbers; (II) Ei(s)≥0, i = 1,2, . . . , n; and (III) On [Sik−1, Sik), Ei(s) is monotone and almost surely differentiable with a positive derivative E0 i(s). To demonstrate it, Figure 1.3.1 shows the effective age for a unit. This process between 0 (or S0) and S2is E(s) = s. At the first event S1, the unit is “minimally repaired”. Thus, in this case patients have no improvement after first occurrence. Repair concept can be translated to biomedical context as “effect of treatment” because we cannot decide if a patient will receive either a “perfect” or a “minimal” repair. Physicians try to do the best for the patient and after the treatment we can observe if the intervention was “perfect” (e.g, patient achieves a complete remission) or “minimal” (e.g., disease is still active). A “perfect repair” occurs to the unit at the second failure time S2. For medical people the treatment is also ”perfect” since the disease disappears, although it is still in the patient but not active. After that, the effective age may be 10 Chapter 1 Introduction: State-of-the-art Figure 1.4: Effective age process, E(s), for an hypothetical case. represented, in general, by a nonlinear function between S2and S3and also between S3and S4. The unit is repaired after the third event between a minimal and perfect repair. In medicine we can say that the patient after third relapse suffers a little improvement produced by the treatment. The fourth failure, S4, is not observed since the observational period for this unit is less than the calendar time. In consequence T4is censored at τ. 1.3.2 Reliability models Next, we describe different existing models which take into account the effect of performed interventions. Most of them are models from reliability or engineering settings. We try to illustrate, using some real examples, how these models can be applied in biomedical or public health problems. Minimal Repair Model In this model, the repair restores the system to its state just before failure. Under this model, the hazard of event occurrence is identical to the intensity just prior to the event occurrence. In this case the effective age corresponds to calendar time, e.g., E(s) = s(in biomedical setting 1.3 Estimation of the survival function 11 this model is called total time model). Thus, in reliability terms one says that the subject is ”minimally repaired” through the intervention. This model has been considered in Brown and Proschan (1983) and Lawless (1987). In biomedical or public health settings we can say that treatment (intervention) is not effective, e.g. patient suffers a minimal or null disease improvement and the disease continues in a stable manner. Some examples can be found involving medical data. For example, let us suppose that patients are involved in a clinical trial. Let us also suppose that we are interested in estimating the time until a headache and that patients are treated with a new drug or placebo after each occurrence. The effective age for patients who received placebo must be calendar time (or total time) since the time until next headache will have the same probability distribution than the last headache occurrence. This is true because patients in the placebo arm are not intervened upon relapse since they do not receive any drug (of course, we omit placebo effect). Perfect Repair Model Upon failure, the failed system is replaced by a new one stochastically identical to the original so that T1, T2, . . . are i.i.d. according to F. These model can be seen as a gap time model where E(s) = s−SN†(s)(where s−SN†(s) represents the elapsed time since the last event occurrence). In medicine one can find some examples when intervention completely treats the disease. A good example is recurrence in some cancers. For example, patients diagnosed to superficial bladder cancer can obtain a complete remission of their disease by means of a surgery intervention after each relapse. It is possible to get a complete remission upon each relapse because we are dealing with an indolent disease which can be controlled easily by means of the surgery. In this case, the perfect repair can be understood as some markers for the disease disappearance. At this moment the disease is not active, but disease is still in the patient since probably disease appears again as a new relapse or recurrence. Brown and Proschan (1983) (BP) generalized the minimal repair model by allowing two types of repairs. After each failure, the system can be perfectly repaired with probability pand, with probability 1 −p, a minimal repair is performed. If the success probability pis made to depend on the time of event occurrence, the Block et al. (1985) (BBS) model is obtained. Hollander et al. (1992) and Presnell, Hollander, and Sethuraman (1994) gives more details of this model. It is difficult to find a biomedical example where the effective age follows this model. 18 Chapter 1 Introduction: State-of-the-art are i.i.d. with common conditional survivor function ¯ F(t|Zi=z) = [ ¯ F0(t)]z= exp µ−zZt 0 λ0(u)du¶, where λ0(·) is the hazard function associated with a baseline survivor function ¯ F0(·). The frailties Z1, Z2, . . .,Znare assumed to be i.i.d. from an unknown distribution function H. In general, the Z’s are not observed, so we are interested in estimating the marginal survivor of Tij, which under this model is given by ¯ F(t) = E{exp(−Z1Λ0(t))}(1.14) where Λ0(t) = −log[ ¯ F0(t)] is the cumulative hazard function of ¯ F0. A common choice of the unknown frailty distribution His a gamma distribution with shape and scale parameters both equal to an unknown parameter αin order to guarantee identifiability. In this case, the common marginal survivor function ¯ Fin (1.14) becomes ¯ F(t) = ·α α+ Λ0(t)¸α .(1.15) The parameter αcontrols the degree of association between interoccurrence times within subject. In particular, as αincreases (decreases), association between interoccurrence times decreases (increases). Letting α−→ ∞, we obtain a model with independent interoccurrence times in which the Tij has a common survivor function of ¯ F0. Pe˜na et al. (2001) showed that the estimation of αand Λ0of (1.15) can be obtained via the maximization of the marginal likelihood function of αand Λ0(·) and with an implementation of the expectation-maximization (EM) algorithm (see, for details, Pe˜na et al., 2001). This estimator of (1.15) is of form ˜ ¯ F(s, t) = "ˆα ˆα+ˆ Λ0(s, t)#ˆα where ˆ Λ0(s, t) is an estimator of the marginal cumulative hazard function Λ0(t). As pointed out in Pe˜na et al. (2001), model parameter estimation can be carried out using the ideas in Nielsen et al. (1992) who suggest that the EM algorithm proposed by Dempster et al. (1977) can be used to obtain the maximizing values (further details in how to apply the EM algorithm when we deal with frailty models may be found in subsection 1.4.2). The shared frailty model can be extended to a model with covariates by means of a multiplicative regression in which frailties act multiplicatively on hazard like in the Cox model. Thus, 1.4 Within-subject correlation 19 proportional hazard model for subject iis written as: λi(s|Zi, Xi) = Ziλ0(s) exp{β0Xi(s)}(1.16) where Ziis the frailty for the ith subject. In applied work, the most widely parametric distribution assumed for frailties is the gamma distribution (Vaupel et al., 1979; Clayton and Cuzick, 1985; Klein, 1992; and Andersen et al., 1993). However, another distributions have been employed (see for instance Hougaard, 1994). Maximum likelihood estimation in the semiparametric frailty model (with gamma-distributed frailties) may be performed using EM algorithm (Dempster et al., 1977) as suggested Gill (1985). This method was further discussed by Nielsen et al. (1992), Klein (1992), and Guo and Rodriguez (1992). Assuming a parametric model, another possibility is to make direct use of the observed data (partial) likelihood (Aalen, 1988). McGilchrist and Aisbett (1991) and McGilchrist (1993) use partial likelihood procedures assuming that frailties follow a log-normal distribution. Penalized likelihood methods has also been studied by several authors: Verwerj and Houwelingen (1994), Therneau et al. (2003), Therneau and Grambsch (2000, section 9.6), Ripatti and Palmgren (2000). Next sections illustrate EM and penalized approaches to make inference for the frailty model. EM algorithm for frailty models We assume that the Z’s are i.i.d. from a distribution H(· | ξ). A common choice of His the gamma distribution with mean 1 and variance 1/ξ, e.g. H= Γ(ξ, ξ). Imposing this restriction is needed to have identifiability as we have pointed out in Subsection 1.4.2. Thus, if Ziare known, the complete log-likelihood is given by n X i=1 µZ∞ 0 Yi(t)[log(λ0(t)) + log(Zi) + β0Xi(t)]dNi(t)−Z∞ 0 Yi(t)Ziexp{β0Xi(t)}λ0(t)dt + log f(Z;θ)¶. As Zcan be viewed as missing data, the problem can be approached using the EM algorithm implemented by Nielsen et al. (1992) in counting process frailty models. The main ingredients of this algorithm are two steps: E (expectation) and M (maximization). We outline briefly both steps, see Andersen et al. (1993) or Nielsen et al. (1992) for further details. Given (Λ0(·), β), the observed times, and the covariates, the conditional expectation of Ziis E-step: E{Zi|Λ0(·), β)}=ξ+Ni(s∗) ξ+Rs∗ 0Yi(v) exp{β0Xi(v)}λ0(v)dv, 20 Chapter 1 Introduction: State-of-the-art where s∗denotes an upper limit of observation times. In the M-step Λ0(·), α, β can be estimated using procedures describe in Appendix A because after the E-step Z’s are set to ˆ Zi. Thus, for example, M-step: ˆ Λ(t) = Zt 0 dN(u) Pn i=1 ˆ ZiYi(u) exp{β0Xi(u)}. To estimate ξAndersen et al. (1993, Subsection IX.4.2) propose to maximize the marginal profile likelihood for ξ. Penalized likelihood estimation The penalized regression formulation for the frailty model is easily developed by making the change Zj= exp(zj). Thus, equation (1.16) converts to λi(s|z, Xi) = λ0(s) exp{β0Xi(s) + z0Mi}, where z0is a vector of frailties and Mis a matrix of nindicator variables such that Mij = 1 when observation iis a reoccurrence of individual jand 0 otherwise. The penalized likelihood method was introduced by Good and Gaskins (1971) in the context of nonparametric probability density estimation. Its use in Cox regression model estimation was proposed by several authors: Zucker and Karr (1990), McGilchrist and Aisbett (1991), Verwerj and Houwelingen (1994), Therneau et al. (2003), and Ripatti and Palmgren (2000). The idea of this method is to maximize penalized partial likelihood equation PPL =l(β, z)−g(z|θ) (1.17) over both βand z. In this equation l(β, z) is the log of the Cox partial likelihood given by l(β, z) = n X i=1 Z∞ 0 Yi(t)(β0Xi(t) + z0Mi)−log   n X j=1 Yj(t) exp{β0Xj(t) + z0Mj}  dNi(t), and gis a penalty function chosen by the investigator to restrict the values of z. The parameter θis a tuning constant which may be pre-specified or adapted to the data to control the amount of shrinkage. Typically, we are interested in choosing the penalty function to “shrink” ztoward zero. Using Newton-Raphson method we can estimate βand zsolving the score equations. The penalty function does not involve β, so we can compute ∂PPL/∂β using ∂PL/∂β which is the 1.4 Within-subject correlation 21 usual partial likelihood for the Cox model. Therefore, the score equation for βare the same as those for Cox model but incorporating zas an offset term (see Therneau and Grambsch, 2000 for further details). In addition to the score vectors, the maximization algorithm requires the Hessian of the penalized partial log-likelihood which is given by H=I+ 0 0 0g00  where I=I(β, z) is the second derivative matrix of partial likelihood, also called information matrix, and g00 denote the second derivative of g. 1.4.3 Cox extension models Many survival models based on Cox proportional hazards have been proposed that handle multiple event data (see Therneau and Grambsch, 2000, Chapter 8; Therneau and Hamilton, 1997; Kelly and Lim, 2000; or Barai and Teoh, 1997 for excellent reviews of these models; Barcel´o, 2002 has a review published in a Spanish journal). In general, there exists two different approaches that extend the Cox model: Wei, Lin, and Weissfeld (1989) (WLW) and Lee, Wei, and Amato (1992) (LWA) marginal, and Prentice, Williams, and Petersen (1981) (PWP) conditional models. On the other hand, frailty models can also be seen as a Cox-based model. Before starting to illustrate both models we mention that there exists another approximation known as the AG model (Andersen and Gill, 1982). The coefficient estimates using AG apporach are exactly the same as those obtained using Cox model, only the standard errors are different. The AG model is the most simple variance corrected model, incorporating robust variance estimators. However, this model requires the strongest assumptions. The main hypothesis of AG model is that repeated events within-subject are independent (given the covariates). This assumption is called “independent increment”, e.g., one event is not affected by previous events. This restriction means that event dependence cannot be estimated with this model, e.g, the model assumes that events do not change the subject and that the subject does not “learn” from previous events. In addition, AG model does not allow one to investigate effects that might change based on event specific covariate effects. However, we have noticed that there exists the possibility to incorporate event dependence via time-dependent covariates. Given these limitations, AG model is recommended when there is no event dependence and no covariate/event effects. The hazard function 22 Chapter 1 Introduction: State-of-the-art of an individual ifor the kevent is given by λik(s|Xik) = λ0(s) exp{β0Xik(s)}. In this case the set of subjects at risk, (e.g. the risk indicator) is given by Yik(s) = I{Si,k−1< s ≤Sik}.(1.18) Marginal models Wei, Lin, and Weissfeld (1989) (WLW) illustrated the marginal model with bladder cancer data set with multiple relapses per patient (see Chapter 2 for a more detailed description of these data). Their method model the marginal distribution of each failure time and no particular structure of dependence among distinct failure times on each subject is imposed. Each recurrence is modelled as a different strata. Data are used in each strata as marginal data, and as Therneau and Hamilton (1997) pointed out, “what would result if the data recorder ignored all information except the given event (type)”. This model is marginal with respect to the risk set since each patient is at risk from the beginning of study and can be at risk for several events simultaneously. The intensity or hazard function for the kth event for the ith subject is λik(s|Xik) = λ0k(s) exp{β0Xik(s)}, For this model, the set of subjects at risk just prior to time sis with respect to the kth event. So, Yik(s) can be given by Yik(s) = I{Sik ≥s}, which corresponds to total time formulation. The estimates of WLW can be either event-specific or overall. The overall estimate proposed by WLW is the weighted average of the event-specific estimates, ˆ β1, . . . , ˆ βk, such that the corresponding weighted average of the robust variance is the smallest possible (Wei et al., 1989). The WLW model presents important disadvantages. Cook and Lawless (1997) pointed out that WLW model is valid only under independent censoring. This disadvantage, however, does not have problems in practice, except if a recurrence is terminal. That is, if the end of study is related to interoccurence times. This model also requires that the data has a maximum number of events. This is a limitation if event-specific estimates become unreliable. 1.4 Within-subject correlation 23 Conditional models The WLW limitations can be solved using a conditional approach. In contrast to AG model, Prentice, Williams, and Petersen (1981) (PWP) propose a conditional model which allow for event dependence via stratification by event number; different events can have different baseline hazards. The main difference with marginal models is that in conditional models a subject cannot be at risk for the n-th event until the (n−1)-th event occurs, hence the name conditional model. Oakes (1991) argues for the conditional approach, and states that the marginal method is inefficient. The conditional model is another variance corrected model and as we have illustrated, it has intuitive appeal because it preserves the order of sequential events in the creation of the risk set and therefore incorporates events dependence. PWP models can be estimated with the data organized in elapsed time (PWP-TT) (i.e., total time risk set or time from each unit’s entry into the observation set) or interocurrence/gap time (PWP-GT) (i.e., gap time risk set or time since the previous event). Thus, the hazard function only differs in the risk intervals formulation. For PWP-TP model the hazard function is given by: λik(s|Xik) = λ0k(s) exp{β0Xik(s)}, and for PWP-GT model by λik(s|Xik) = λ0k(s−sk−1) exp{β0Xik(s)}. In consequence, Yik(s) is different for each formulation. In the PWP-TT model Yik(s) corresponds to Yik(s) = I{Sik ≥s}, and for PWP-GT model it is given by Yik(s) = I{Tik > s}. The choice between PWP-TT or PWP-GT depends on whether we are interested in the time that has elapsed since the patient entered the study or since the last recurrence. We notice that PWP model is a stratified AG model. Some of these previous models have been compared using real and simulated data, giving different results as it is illustrated in Wei and Glidden (1997), Gao and Zhou (1997), Clayton (1994), Lin (1994), Therneau and Hamilton (1997), Barai and Teoh (1997), and Therneau and Grambsch (2000). The AG model is maybe the most used because of its efficiency as Therneau and Grambsch (2000, pag 229) conclude from their hidden covariate Monte Carlo simulation. In 24 Chapter 1 Introduction: State-of-the-art addition, in the case of variance corrected models, the AG model gives the most reliable estimates of the overall effect. However, as Kelly and Lim (2000) pointed out, it remains unclear which models are suitable for recurrent event data, as well as the differences between existing models. They do not recommend LWA model because it allows a subject to be at risk several times for the same event. The WLW model overestimates treatment effect and they do not recommend it. Finally, the authors propose to use PWP-GT for analyzing recurrent event data. However, when there exists within-subject correlation, they recommend to use methods different from those which are based on a robust variance estimation, like random effects models. Software All Cox extension models, outlined previously, can be fitted using both S-plus (MathSoft, 1997) and R (R Development Core Team, 2005, Ihaka and Gentleman, 1996) functions creating an appropriate data set. After obtaining these data sets, we can fit the models mentioned above using coxph function to fit Cox models and the functions strata,cluster, or frailty. Therneau and Grambsch (2000, Chapter 8) or Therneau and Hamilton (1997) explain how to fit any of the models described in the previous sections. An excellent review of the different software packages used for analyzing correlated survival data has been recently written by Kelly (2004). 1.5 General class of models As we have previously illustrated, there are currently several models and methods of analysis used for recurrent event data (see for instance Hougaard, 2000; or Therneau and Grambsch, 2000, Chapter 8). However, as Pe˜na and Hollander (2004) point out, there is still a need for a general and flexible class of models that simultaneously incorporates the effects of covariates or concomitant variables, the impact on the unit of accumulating event occurrences, the effect of latent or unobserved variables which, for each unit, endow correlation among the inter-event times, as well as the effect of performed interventions after each event occurrence. Most existing extensions of the Cox model deal with the majority of these effects. However, as we have pointed out in Section 1.3.1 these models assume that the effect after each intervention is always the same (minimal or perfect intervention if risk interval is total or gap time, respectively). As we have also illustrated, in many biomedical settings not always the effect of intervention is 1.5 General class of models 25 the same. So, we need to use a more general model that allow us to incorporate different effects that affect to event occurrences. Pe˜na and Hollander (2004) proposed a new class of models which generalize most of existing reliability and Cox-based models. This new model will allow us to connect reliability models, which models the effect of intervention via the effective age, with biomedical models, which incorporates the effect of concomitant covariates and the correlation among interoccurrence times. This model is briefly described in the next section. 1.5.1 Pe˜na and Hollander model Let Z= (Z1, Z2, . . . , Zn) be a vector of independent and identically distributed (i.i.d.) positivevalued random variables from a parametric distribution H(z;ξ) = Pr(Z≤z|ξ) where ξis a finitedimensional parameter taking values in Ξ ⊆ <r. These variables are unobservable random factors affecting the event occurrences for the subjects. Also, let F={Fs: 0 ≤s≤s∗}be a filtration or history on some probability space (Ω,F,P) such that Xis and Y† is are predictable and such that N† is are counting processes with respect to F. Finally, let {Ei(s) : 0 ≤s≤s∗}, i = 1,2, . . . , n, the effective age processes satisfying the conditions described in Section 1.3.1. The class of models is obtained as follows. Conditionally on Z, the F-compensator of N† iis {A† i(s|Z,Xi) : 0 ≤s≤s∗}with A† i(s|Z,Xi) = Zs 0 Y† i(v)λi(v|Z,Xi)dv, (1.19) where λi(s|Z,Xi) = Ziλ0[Ei(s)] ρ[N† i(s−); α]ψ[β0Xi(s)].(1.20) This means that the process M† i(s|Z,Xi) = N† i(s)−A† i(s|Z,Xi) is a square-integrable F-martingale. In (1.20), λ0(·) is an unknown baseline hazard rate function. The effects of accumulating event occurrences are encoded in ρ(·;α) : Z+≡ {0,1,2, . . .}→<+which has a known functional form with ρ(0; α) = 1 and with α∈ A ⊆ <p. The effect of covariates are considered in ψ(·) which is a nonnegative link function of known functional form with β∈ B ⊆ <q. The dependence between interoccurrence times is modelled with Ziwhich are unobserved frailties. The model also incorporates the effect of performed interventions, Ei(s), via the baseline hazard function. Thus, the unknown model parameters are (λ0(·), α, β, ξ), where λ0(·) is non-parametrically specified, and α,β, and ξare finite-dimensional parameters. 26 Chapter 1 Introduction: State-of-the-art The statistical identifiability of this class of models without frailties has been established in Theorem 1 of Pe˜na and Hollander (2004). The authors also showed that this class of models subsumes many existing models in the literature (see Pe˜na and Hollander, 2004). In particular, some of models used in biomedical settings are special cases of this general model as we will illustrate in Chapter 5. 1.6 Thesis Overview In the following chapters, we present the work that we have been developing during last few years. In the Chapter 2 we describe four data sets analyzed for the elaboration of this dissertation. We would like to emphasize that two of them were obtained from the Institution where I was working until July 2005. These data sets were created specially for this PhD thesis. The other two come from medical literature and have been widely used when researchers propose new methods to analyze recurrent event data. I must say that my motivation was to find real examples, which were suitable to be analyzed with the statistical methods I was working on. Perhaps the most relevant analysis was the one about lymphoma relapses, in which we included the information about the treatment after each relapse by means of the effective age, which had never been done before. Then, in Chapter 3, we develop some procedures for estimating confidence intervals for median survival time or, in general, for some quantile. There, we propose some asymptotic confidence intervals which are based on asymptotic variances from existing estimators for survival function when we deal with recurrent events. We also propose how to estimate these confidence intervals using bootstrap techniques. The main contribution of this chapter has been the examination of the question of how to do bootstrapping in the presence of recurrent event data arising from a sum-quota data accrual scheme and informativeness of right-censoring mechanism. We show how to get bootstrap samples from the observed data, as many people normally do, is not correct when data are correlated. Chapter 4 deals with procedures for estimating the parameters for the general model for recurrent events proposed by Pe˜na and Hollander (2004). One possibility is to use the EM algorithm as in a joint work with Professors Edsel A. Pe˜na and Elizabeth E. Slate showed (paper currently in first revision in Journal of Statistical Planning and Inference and published as a 1.6 Thesis Overview 27 technical report in Pe˜na, Slate, and Gonz´alez, 2003). This work is presented in the Appendix B since my contribution was the programming of all the procedures exposed, as well as proposing alternative maximization methods. A result of that is the gcmrec package and this is why it is include in the thesis. These functions allowed us to perform simulation studies, and analyze the examples shown on the publication where I also participated. We have included this work in the appendix of the thesis because Chapter 4 is based on the notation and the results we get in it. It is well-known that the EM algorith method have general drawbacks such as neither estimates of the variance of parameter nor frailty are directly estimated. Thus, still in Chapter 4, we proposed two alternative approaches, based on penalizing the likelihood, to fit Pe˜na and Hollander model. One of them follows Therneau et al.’s 2003 work. They proposed to penalize regression coefficients. This approach still continues to have problems because the convergence can be slow and the variance of frailty cannot be directly estimated. Then, we propose to adopt another method of penalization described in Rondeau et al. (2003). Their idea is to penalize the full likelihood, instead of the partial likelihood as Therneau et al. (2003) proposed, and to obtain smooth estimates of the hazard function. This method has the main advantage of giving an estimate of variance of the frailty variance. Finally, Chapter 5 address the problem of how to incorporate the effective age process in biomedical settings. So far, this concept has only be used in reliability problems. Our main contribution in this chapter has been to illustrate how to use the information regarding the effects of treatments or interventions after cancer relapses for modelling the effective age. Our motivation was firstly due to the fact that by analyzing some data sets, and carrying-out some simulations, we showed that one can obtain different results by using different effective age formulations. In addition, as some physicians pointed out“it is necessary a model designed specifically for relapsing patients” (MacLaughlin, 2002). A section reproducing the on-line documentation of Rpackages developed for this dissertation that also can be obtained at http://www.r-project.org/ are described in the Appendix D. We end the chapters by showing how to analyze the data presented in each chapter using these R functions. 34 Chapter 2 Studies with Recurrent Event Data 2.2 Non-Hodgkin’s Lymphoma Cancer Relapses The indolent non-Hodgkin’s lymphomas (NHL) constitute a heterogeneous group of lymphoproliferative disorders. They encompass what were called low grade and some categories of intermediate grade NHL in the Working Formulation (Cheson et al., 1999). They are categorized based on pathologic and cytologic features. The indolent lymphomas include different subtypes of lymphomas such as Follicular Lymphomas, Small Lymphocytic Lymphoma, Lymphoma Marginal Zone, or Sezary Syndrome among others. Low grade lymphomas are associated with relatively prolonged survival. Because it is considered an indolent, but not curable, type of cancer, patients tend to relapse over time. Thus, patients are treated after each recurrence with intensive therapeutic approaches in an attempt to increase the time until next relapse (that is, to increase the disease-free survival). The treatments may produce different responses (e.g., complete response, CR, partial response PR, or null response NR) depending on disease status after therapy. It is well known that these responses may modify the probability of a subsequent relapse, and hence this intervention effect should be taken into account when modelling this type of data. The data consist of the times to relapse, in months, for 63 patients with clinical, histopathological, and immunophenotypes of primary cutaneous marginal zone B-cell lymphoma (PCMZCL) as a particular subtype of indolent lymphoma. An analysis of a subset of these data based on 22 patients with a specific subtype of cutaneous lymphoma was presented in a recent paper Servitje et al. (2002). We use the date of first treatment as the beginning of the study. Relapsing times were considered as the difference between last relapse and the current one. Figure 2.2 shows a graphical representation of recurrences for this data set. 2.2.1 Summary of the data set We have also obtained information about the response achieved after treatment upon relapses (CR, PR, or NR), depending on the disease status, for each relapse for each of the 63 patients. The total number of relapses among all patients is 112. The fraction of patients with no relapse is 57%, and only 7% have 3 or more events. The median follow-up time is 2.9 years (range 1 month to 13.5 years). Thirty eight of 49 (77.8%) responses to treatments administered after relapses are CRs, 9 (14.3%) are PRs, and 2 (4.5%) NRs. This information will be useful to model the effective age as we will illustrate along this thesis. As we can see in the Figure 2.3 the type of response to 2.2 Non-Hodgkin’s Lymphoma Cancer Relapses 35 0 50 100 150 0 10 20 30 40 50 60 Calendar Time Subject Figure 2.2: Graphical representation of the lymphoma data set. The graphic shows the times (◦) and censoring (×) of PCMZCL cancer recurrence for 63 patients. 0 50 100 150 5 10 15 20 Calendar Time Subject Figure 2.3: Graphical representation of the lymphoma data set including information about the response to treatment after relapses. The graphic shows the times (◦: complete remission, M: partial response, ¤: no response) and censoring (×) of non-Hodgkin’s Lymphoma cancer recurrence only for 45 selected patients for for improving the presentation of the graphic. 36 Chapter 2 Studies with Recurrent Event Data the treatment is related to the time to next relapse. We will also include in the analyses the covariates X1: gender of patient (0=Male, 1=Female); X2: delay between first symptom and date of first treatment as a continuous variable (in years); and X3: lesions involved at diagnosis (0=Single, 1=Localized, 2=More than one nodal site, 3=Generalized), encoded as three indicator variables. The distribution of these covariates is as follows: 73 patients are males (65.2%) and 39 females (34.8%); the median time of the delay between first symptom and first treatment is 29.7 months (range 1 to 144 months); 28 patients (25.0%) presented single lesions at diagnosis, 43 localized lesions (38.4%), 35 more than one nodal site (31.2%), and 6 (5.4%) patients had generalized lesions. 2.3 Bladder Cancer Relapses Last data set related to cancer is on recurrences of bladder tumor. These data have been used by many people to demonstrate methodology for recurrent event modelling and they can can be obtained from the survival package (Lumley and Therneau, 2003) in the RLibrary. Wei, Lin, and Weissfeld (1989) analyzed these data using marginal approach. These data provide the times to recurrence of bladder cancer for n= 85 subjects with superficial bladder tumors, which were removed when they entered the study. Forty seven of these patients were randomized into the placebo group, and 38 into the thiotepa group. Many patients have multiple recurrences of tumors in the study, and new tumors were removed at each visit. The data set contains the first four recurrences of the tumor for each patient, and each recurrence time was measured from the patient’s entry time into the study. The total number of recurrences was 112. The covariates are X1, the treatment indicator (1 = placebo; 2 = thiotepa); X2, the size (in cm) of the largest initial tumor; and X3, the number of initial tumors. Figure 2.4 shows a graphical representation of recurrences for this data set. No information about effective age is available for these data. 2.4 MMC data set This data set belong to data from a study concerning small bowel motility (Husebye et al. (1990)). The aim of their analysis is to estimate the mean length of the Migratory Motor Complex (MMC) period (i.e., the mean interoccurrence time between two contractions in the small bowel during the 2.4 MMC data set 37 0 10 20 30 40 50 60 0 20 40 60 80 Calendar Time Subject Figure 2.4: Graphical representation of the bladder data set used by Wei et al. (1989). The graphic shows the times (×) and censoring (◦) of bladder cancer recurrence for 85 subjects. The 38 top subjects are treated with thiotepa and the upper 47 subjects are from placebo group. interdigestive cycle). This data set was analyzed in Aalen and Husebye (1991) using a variance component model and an intensity-based formulation with a gamma frailty component using a parametric Weibull model. Then Pe˜na et al. (2001) also analyze this data using two new estimators which are described in Sections 1.3.3 and 1.4.2 respectively. Aalen and Husebye (1991) stated that“the consecutive MMC periods for each individual appear (to be) approximate renewal process” but we need to verify this assumption. To do so, Pe˜na et al. (2001) suggested that since formal statistical methods for checking this i.i.d. assumption are not yet available, a graphical method may be employed. We will illustrate it in Chapter 3. Chapter 3 Confidence Intervals for Median Survival This chapter addresses the problem about how to construct confidence intervals for the median survival time of a recurrent event. Two different approaches have been employed. One of them is based on asymptotic variance of Pe˜na et al. (2001) and Wang and Chang (1999) estimators (see Sections 1.3.3 and 1.3.4 respectively) and some transformations. The other one uses bootstrap techniques. Two types of recurrent event models are considered: first is a model where the interevent times are independent and identically distributed (see Section 1.3.3), and second is a model where the inter-event times are associated, with the association arising from a gamma frailty model (Section 1.4.2). Both bootstrap and asymptotic confidence intervals are studied through simulation. Weak convergence is proved and asymptotic confidence intervals are built according to these results. On the other hand, one of the major goals of this chapter is to study bootstrapping schemes for estimating the sampling distribution of estimators of the median survival time distribution in the presence of recurrent event data. Another important goal is to determine how to compute pointwise confidence intervals for median survival time using asymptotic theory for PSH and WC estimators. We now discuss the motivations for the present Chapter. Pe˜na et al. (2001) (under i.i.d. model) and Wang and Chang (1999) give asymptotic properties of their estimators. This fact allows us to approximate analytically the variance of survival function. Thus, we can use their estimators to compute confidence intervals for some quantile, ζp, in particular for median survival 39 40 Chapter 3 Confidence Intervals for Median Survival time, ζ1/2. However, the estimator proposed by Pe˜na et al. (2001) under correlated case does not have a closed form for its variance. If we can assume that our data have been generate under a gamma frailty model, we need an alternative method to estimate the variance of median survival time. As we have previously mentioned, we can use resampling techniques. Bootstrap procedures allow us to construct some confidence intervals for bootstrapped percentiles of ζ∗ p. However, before computing this variance, we need to determine how to obtain bootstrap samples by taking into account the true nature of this kind of data. We need to recognize that the informativeness of right-censoring mechanism, and the impact of the sum-quota accrual scheme, may influence the estimation of survival function, and consequently must be considered in generating bootstrap samples. Section 3.1 provides some procedures to obtain asymptotic confidence intervals for median survival time (and other quantiles) using asymptotic results for PSH and WC estimators. Different bootstrapping schemes are described in section 3.2 for the i.i.d. model and for the correlated interoccurrence times model. In Section 3.3 a simulation is used to compare and discuss the statistical properties of median survival time estimated using these different bootstrap plans. In this section, we also compare the results obtained using improved bootstrap plans with those obtained using asymptotic theory. In Section 3.4 we apply these procedures to three real data sets described in Chapter 2. We illustrate the problems that one can have when analyze these kind of data. The first example belongs to data from small bowel motility study. This is a good example when interoccurrence times follow an i.i.d. model. The second example involves hospital readmissions in patients diagnosed with colorectal cancer. This example is an excellent example when interoccurrence times are correlated. The third example analyzes tumor recurrences in bladder cancer. This example shows how the results can be different depending on the model selection. Finally, Section 3.5 shows how to analyze data with recurrent events using survrec package. 3.1 Estimation of median survival time and other quantiles After estimating the survival function of the interoccurrence times, we are usually interested in estimating some quantile of this function, ζp(x). For example, in biomedical settings, one may want to estimate the median survival interoccurrence time, ζ1/2(x). After that, the main 3.1 Estimation of median survival time and other quantiles 41 concern is assessing the variability of the estimator, ˆ ζp, by estimating its standard deviation, and constructing confidence intervals for ζpbased on ˆ ζp. Different approaches can be used to do so. Firstly, we can use asymptotic methods (see Burr and Doss, 1993 or Dabrowska and Doksum, 1987) taking into account that it is more difficult to apply this methodology in the case of ζp(x) than in other situations. Secondly, as several authors have stated, we can estimate the sampling variability of ˆ ζpmore accurately using resampling techniques (see for instance Efron, 1982; Efron, 1985b; Efron, 1985a; Bickel and Frieedman, 1981; Beran, 1982; and Singh, 1981). In particular, some authors have shown that these techniques are useful in the case of survival analysis under the Cox’s model (see Burr, 1994, Hjort; 1985). Let 0 < p < 1, define ζp=F−1(p) = inf{t:F(t)≥p}= inf{t:¯ F(t)≤1−p}(3.1) as the pth quantile of the interoccurrence distribution function, F. The quantiles can be estimated by taking the right-continuous inverse of the nonparametric estimated survival function ˆ ¯ For ˆ S. Let ˆ ¯ F(n)be the nonparametric estimation of the interoccurrence times survival function using either ˆ ¯ F, or ˆ S, estimators defined in (1.8) and (1.11) respectively. Thus, ˆ ζp= inf{t:ˆ ¯ F(n)(t)≤1−p}.(3.2) Proposition 1 Let ζpand ˆ ζpbe defined in (3.1) and (3.2) respectively and let us assume the conditions stated in Theorem 1 from Pe˜na et al. (2001) and in Theorem 1 from Wang and Chang (1999). For the PSH estimator, let also assume the condition y(∞, t∗)>0, for any t∗> ζpwhere y(∞, t∗) = lims→∞ E{Y(s, t)}. Then √n(ˆ ζp−ζp)d −→Nµ0,p2σ2(ζp) f2(ζp)¶ where f=−¯ F0, and σ2corresponds to σ2 PSH or σ2 WC (defined in equations 1.9 and 1.12 respectively) depending on whether we estimate ˆ ¯ F(n)using either ˆ ¯ For ˆ S(respectively). Proof: We demonstrate this proposition using the delta-method illustrated in Andersen et al. (1993, section IV.3.4). We write ζp=φ(¯ F) and ˆ ζp=φ(ˆ ¯ F(n)), where φis the function defined by φ(G) = G−1(p) = inf{x:G(x)≤1−p}(3.3) 42 Chapter 3 Confidence Intervals for Median Survival Using Proposition II.8.4 from Andersen et al. (1993) and f(ζp)>0, we have that φis (tangentially) compactly differentiable at ¯ Fwith the derivative given by ¡dφ(¯ F)¢(p)·h=h(¯ F−1(p)) f(¯ F−1(p)),(3.4) where his continuous at ¯ F−1(p). Using functional delta-method we can approximate √n(φ(ˆ ¯ F(n))−φ(¯ F)) by dφ(¯ F)√n(ˆ ¯ F(n)− ¯ F), where dφ(¯ F) is the derivative of φat ¯ Fand it acts on √n(ˆ ¯ F(n)(t)−¯ F(t)) in a linear way. Now, using this approximation, and the Theorem IV.3.2 of Andersen et al. (1993), if we demonstrate that as n−→ ∞ √n(ˆ ¯ F(n)−¯ F)d −→Z(3.5) where Z=¯ FU, with Ubeing a Gaussian martingale with U(0) = 0 and cov(U(s1), U(s2)) = σ2(s1∧s2), we will be able to conclude that √n(φ(ˆ ¯ F(n))−φ(¯ F)) d −→dφ(¯ F)·Z and then, √n(ˆ ζp−ζp)d −→ pU(¯ F−1(p)) f(¯ F−1(p)) ∼Nµ0,p2σ2(ζp) f2(ζp)¶ Thus, to demonstrate (3.5) we need to prove that (i) √n(ˆ ¯ F−¯ F)d −→Z1 (ii) √n(ˆ S−¯ F)d −→Z2 where Z1and Z2must be a Gaussian martingale with U(0) = 0 and cov(U(s1), U(s2)) = σ2(s1∧s2) as we have previously mentioned. The point (ii) is proved in Theorem 1 from Wang and Chang (1999) while point (i) may be proved using Theorem 2b from Pe˜na et al. (2001) as follows. This Theorem states that for all fixed s∈[0,∞] √n(ˆ ¯ F(s, t)−¯ F(t)) d −→Z. Taking s=∞we have ˆ ¯ F(∞, t)≡ˆ ¯ F(t), and adding the condition y(∞, t∗)>0, where y(∞, t∗) = lims→∞ E{Y(s, t∗)}, √n(ˆ ¯ F(t)−¯ F(t)) d −→Z. (3.6) ¤ 3.1 Estimation of median survival time and other quantiles 43 Let us notice that the condition that y(∞, t∗)>0 limits the interval in which we have convergence to be [0, t∗] (see conditions of Theorem 1 from Pe˜na et al., 2001). Let us also note that it will not work out if t∗=∞. However this is an assumable situation at least in biomedical settings. Using this Proposition, it is easy to estimate the asymptotic variance of the quantile inserting ˆσ2 PSH (1.7) or ˆσ2 WC (1.13) for σ2(ζp) and ˆ ζp(3.2) for ζpin the expression of Proposition 1. In order to estimate the density f=−¯ F0we can use the kernel function estimator. Using an uniform kernel function, we can estimate fby ˆ f(t) = 1 2b(¯ F(n)(t−b)−¯ F(n)(t+b))) The problem, here, is to give some value for b(bandwidth parameter). A conservative choice is to take bto have 50% of observed times in the interval (t−b, t +b). This criteria will be used in the examples. It is well known that bandwidth selection is a difficult problem. For this reason, it is better to apply Brookmeyer and Crowley’s 1982 procedure that do not need to estimate f(see Brookmeyer and Crowley, 1982). Thus, we take as an approximate 100(1-κ)% confidence interval for ζpall values ζ0 pwhich satisfies |g(¯ F(n)(ζ0 p)) −g(1 −p)| |g0(¯ F(n)(ζ0 p)) |g(¯ F(n)(ζ0 p))ˆσ(ζ0 p)≤cκ/2,(3.7) where cκ/2is the upper κ/2 quantile of the standard normal distribution. Brookmeyer and Crowley (1982) considered g(x) = x, but we can also consider some other transformations like g(x) = log(−log(x)) or g(x) = arcsin √n(see for instance, Kalbfleisch and Prentice, 1980 or Thomas and Grunkemeier, 1975). We may also estimate a confidence interval for ζpas follows: Step 0. For a given p, estimate the pth quantile, ˆ ζpusing (3.2). Step 1. Estimate the confidence interval of the survival function, ˆ ¯ F(n), at t=ˆ ζpas follows: 1.1. If g(x) = x, take ˆ ¯ F(n)(t)±cκ/2ˆσ(t)ˆ ¯ F(n)(t) (3.8) 1.2. If g(x) = log(−log(x)), take ˆ ¯ F(n)(t)exp{±cκ/2ˆσ(t)/log(ˆ ¯ F(n)(t))}(3.9) 50 Chapter 3 Confidence Intervals for Median Survival 95% BPCI MSE (bias) EC ˆµ˜µˆσ2 n=15 Plan I 283 (8.1) 88.4 0.19 0.17 1,034 Plan II 291 (11.3) 94.3 0.21 0.19 1,024 Plan III 291 (11.9) 95.2 0.22 0.20 1,119 Plan IV 703 (10.3) 93.6 0.32 0.28 3,291 Plan V 687 (10.4) 94.2 0.32 0.28 3,421 n=50 Plan I 66 (3.1) 93.3 0.10 0.10 75 Plan II 66 (3.7) 94.6 0.10 0.10 63 Plan III 66 (3.7) 94.8 0.10 0.10 65 Plan IV 141 (3.2) 94.9 0.15 0.15 202 Plan V 142 (3.2) 94.5 0.15 0.15 198 n=80 Plan I 39(2.1) 94.1 0.08 0.08 34 Plan II 38(2.5) 95.4 0.08 0.08 28 Plan III 38(2.6) 95.2 0.08 0.08 28 Plan IV 84(1.9) 95.3 0.12 0.12 88 Plan V 84(1.9) 95.4 0.12 0.12 91 Table 3.1: Simulation results for 2,000 samples and 500 bootstrap replicates under the i.i.d. model. Mean square error (MSE) (×105) and proportion of MSE due to bias. Empirical Coverage (EC) and mean (ˆµ), median (˜µ) and variance (ˆσ2) of the length (×105) of 95% bootstrap percentile confidence intervals (BPCI). Results for the first five bootstrap schemes, varying sample sizes, θ= 1/3 and ν=1. 3.3 Simulation Study 51 95% BPCI MSE (bias) EC ˆµ˜µˆσ2 n=15 Plan I 140 (6.2) 90.4 0.16 0.15 926 Plan II 155 (9.8) 95.3 0.19 0.17 915 Plan III 156 (9.9) 95.1 0.20 0.18 965 Plan IV 321 (8.9) 94.6 0.25 0.24 1,562 Plan V 335 (9.1) 95.1 0.25 0.24 1,635 n=50 Plan I 42 (2.8) 94.3 0.07 0.07 52 Plan II 43 (3.0) 95.6 0.07 0.07 50 Plan III 40 (3.0) 94.9 0.07 0.07 52 Plan IV 111 (2.9) 94.3 0.11 0.11 128 Plan V 115 (3.0) 94.7 0.11 0.11 135 n=80 Plan I 25(2.0) 94.5 0.05 0.05 26 Plan II 26(2.1) 95.4 0.05 0.05 20 Plan III 28(2.2) 95.6 0.05 0.05 21 Plan IV 44(2.0) 95.1 0.10 0.10 56 Plan V 45(2.0) 94.8 0.10 0.10 59 Table 3.2: Simulation results for 2,000 samples and 500 bootstrap replicates under the i.i.d. model. Mean square error (MSE) (×105) and proportion of MSE due to bias. Empirical Coverage (EC) and mean (ˆµ), median (˜µ) and variance (ˆσ2) of the length (×105) of 95% bootstrap percentile confidence intervals (BPCI). Results for the first five bootstrap schemes, varying sample sizes, θ= 1/6 and ν=1. 52 Chapter 3 Confidence Intervals for Median Survival 95% BPCI MSE (bias) EC ˆµ˜µˆσ2 α=2 n=15 Plan IV 2,576 (15.2) 93.4 0.57 0.45 16,985 Plan V 2,362 (15.6) 93.8 0.58 0.46 16,935 Plan VI 1,876 (7.2) 92.1 0.45 0.34 13,115 Plan VII 1,885 (7.3) 92.3 0.45 0.34 13,506 n=50 Plan IV 456 (5.4) 94.2 0.26 0.24 15,422 Plan V 456 (5.5) 94.2 0.26 0.23 16,195 Plan VI 258 (0.1) 93.8 0.20 0.19 518 Plan VII 256 (0.1) 94.1 0.20 0.19 526 n=80 Plan IV 265 (5.6) 94.8 0.20 0.19 426 Plan V 267 (5.8) 95.1 0.20 0.19 434 Plan VI 167 (0.2) 94.4 0.16 0.15 206 Plan VII 168 (0.2) 94.8 0.16 0.15 200 α=6 n=15 Plan IV 1,203 (12.7) 92.9 0.40 0.33 7,160 Plan V 1,154 (12.9) 93.1 0.41 0.33 7,085 Plan VI 927 (16.7) 92.2 0.31 0.26 3,831 Plan VII 921 (16.9) 92.7 0.31 0.26 4,667 n=50 Plan IV 209 (4.0) 93.7 0.18 0.17 319 Plan V 208 (4.1) 93.9 0.18 0.17 352 Plan VI 127 (2.0) 93.7 0.14 0.13 153 Plan VII 127 (2.0) 93.5 0.14 0.13 151 n=80 Plan IV 120 (3.7) 95.3 0.14 0.14 157 Plan V 121 (3.9) 95.4 0.14 0.14 158 Plan VI 72 (1.0) 95.3 0.11 0.10 78 Plan VII 72 (1.0) 95.2 0.11 0.10 77 Table 3.3: Simulation results for 2,000 samples and 500 bootstrap replicates under a gamma frailty model with shape and scale parameter set equal to α=2 and α=6. Mean square error (MSE) (×105) and proportion of MSE due to bias. Empirical Coverage (EC) and mean (ˆµ), median (˜µ), and variance (ˆσ2) of the length (×105) of 95% bootstrap percentile confidence intervals (BPCI). Results for the last four bootstrap plans, varying sample sizes, θ= 1/3 and ν=1. 3.3 Simulation Study 53 95% BPCI MSE (bias) EC ˆµ˜µˆσ2 α=2 n=15 Plan IV 1,562 (12.2) 94.4 0.30 0.20 8,965 Plan V 1,234 (12.1) 93.8 0.32 0.21 8,435 Plan VI 976 (6.4) 93.1 0.25 0.24 7,215 Plan VII 985 (6.3) 93.3 0.25 0.24 7,560 n=50 Plan IV 206 (4.6) 94.4 0.14 0.14 7,922 Plan V 207 (4.2) 94.5 0.24 0.13 7,690 Plan VI 152 (1.4) 94.6 0.10 0.09 218 Plan VII 153 (1.1) 94.5 0.10 0.09 261 n=80 Plan IV 151 (4.7) 94.5 0.10 0.09 261 Plan V 157 (5.0) 94.8 0.10 0.09 243 Plan VI 107 (1.2) 95.2 0.08 0.07 106 Plan VII 108 (0.9) 95.1 0.08 0.07 101 α=6 n=15 Plan IV 842 (9.7) 93.9 0.22 0.19 3,262 Plan V 758 (9.5) 93.8 0.23 0.20 3,381 Plan VI 469 (10.2) 93.1 0.16 0.11 1,831 Plan VII 450 (10.8) 93.4 0.17 0.11 2,005 n=50 Plan IV 112 (3.1) 94.7 0.08 0.07 157 Plan V 111 (3.2) 94.2 0.08 0.07 145 Plan VI 75 (1.6) 94.5 0.06 0.06 79 Plan VII 77 (1.5) 94.0 0.06 0.05 65 n=80 Plan IV 70 (3.1) 95.1 0.06 0.06 89 Plan V 69 (2.9) 95.2 0.06 0.06 82 Plan VI 38 (1.0) 94.8 0.04 0.04 35 Plan VII 38 (1.1) 94.3 0.04 0.04 26 Table 3.4: Simulation results for 2,000 samples and 500 bootstrap replicates under a gamma frailty model with shape and scale parameter set equal to α=2 and α=6. Mean square error (MSE) (×105) and proportion of MSE due to bias. Empirical Coverage (EC) and mean (ˆµ), median (˜µ), and variance (ˆσ2) of the length (×105) of 95% bootstrap percentile confidence intervals (BPCI). Results for the last four bootstrap plans, varying sample sizes, θ= 1/6 and ν=1. 54 Chapter 3 Confidence Intervals for Median Survival PSH estimator WC estimator EC ˆµ˜µˆσ2EC ˆµ˜µˆσ2 n=15 Estimating f(t)a93.6 0.19 0.17 479 88.5 0.27 0.23 2,458 B&C no-trans. 99.9 0.40 0.38 1,685 99.8 0.63 0.58 6,922 B&C log-log-trans. 94.9 0.18 0.17 427 90.5 0.26 0.24 1,282 B&C arcsin-trans. 94.8 0.18 0.17 408 90.8 0.25 0.23 1,255 Bootstrap (Method 2) 94.3 0.21 0.19 1,024 93.6 0.32 0.28 3,291 n=50 Estimating f(t) 93.7 0.10 0.10 28 94.2 0.15 0.14 123 B&C no-trans. 99.9 0.20 0.20 146 99.9 0.32 0.31 574 B&C log-log-trans. 94.6 0.10 0.10 54 94.0 0.15 0.14 186 B&C arcsin-trans. 94.4 0.10 0.10 53 94.3 0.15 0.14 183 Bootstrap (Method 2) 94.6 0.10 0.10 63 94.9 0.15 0.15 202 n=80 Estimating f(t) 94.2 0.08 0.08 11 93.2 0.11 0.11 48 B&C no-trans. 99.9 0.16 0.16 65 99.9 0.24 0.24 245 B&C log-log-trans. 95.1 0.08 0.08 25 94.0 0.12 0.11 81 B&C arcsin-trans. 94.8 0.08 0.08 25 94.0 0.12 0.11 79 Bootstrap (Method 2) 95.4 0.08 0.08 28 95.3 0.12 0.12 88 aestimating f(t) in the expression of asymptotic variance given in Proposition 1 Table 3.5: Simulation results for 2,000 samples under the i.i.d. model. Empirical Coverage (EC) and mean (ˆµ), median (˜µ) and variance (ˆσ2) of the length (×105) of 95% different pointwise Brookmeyer and Crowley’s (B&C) confidence intervals based on asymptotic variance of both PSH and WC estimators. Results for selected sample sizes, θ= 1/3 and ν=1. 3.3 Simulation Study 55 PSH estimator WC estimator EC ˆµ˜µˆσ2EC ˆµ˜µˆσ2 n=15 Estimating f(t)a93.6 0.19 0.17 479 88.5 0.27 0.23 2,458 B&C no-trans. 99.9 0.40 0.38 1,685 99.8 0.63 0.58 6,922 B&C log-log-trans. 94.9 0.18 0.17 427 90.5 0.26 0.24 1,282 B&C arcsin-trans. 94.8 0.18 0.17 408 90.8 0.25 0.23 1,255 Bootstrap (Method 2) 95.3 0.19 0.17 915 94.6 0.25 0.24 1,562 n=50 Estimating f(t) 93.7 0.10 0.10 28 94.2 0.15 0.14 123 B&C no-trans. 99.9 0.20 0.20 146 99.9 0.32 0.31 574 B&C log-log-trans. 94.6 0.10 0.10 54 94.0 0.15 0.14 186 B&C arcsin-trans. 94.4 0.10 0.10 53 94.3 0.15 0.14 183 Bootstrap (Method 2) 95.3 0.19 0.17 63 94.3 0.11 0.11 128 n=80 Estimating f(t) 94.2 0.08 0.08 11 93.2 0.11 0.11 48 B&C no-trans. 99.9 0.16 0.16 65 99.9 0.24 0.24 245 B&C log-log-trans. 95.1 0.08 0.08 25 94.0 0.12 0.11 81 B&C arcsin-trans. 94.8 0.08 0.08 25 94.0 0.12 0.11 79 Bootstrap (Method 2) 95.4 0.08 0.08 28 95.1 0.10 0.10 56 aestimating f(t) in the expression of asymptotic variance given in Proposition 1 Table 3.6: Simulation results for 2,000 samples under the i.i.d. model. Empirical Coverage (EC) and mean (ˆµ), median (˜µ) and variance (ˆσ2) of the length (×105) of 95% different pointwise Brookmeyer and Crowley’s (B&C) confidence intervals based on asymptotic variance of both PSH and WC estimators. Results for selected sample sizes, θ= 1/6 and ν=1. 56 Chapter 3 Confidence Intervals for Median Survival 3.4 Examples In this section we will compute the confidence intervals for median survival time and its confidence interval using MMC, readmission and bladder data sets. We illustrate how to estimate PSH and WC using the survrec package for fitting these data in the Section 3.5. 3.4.1 MMC data set The first data set pertains to data from the study concerning small bowel motility studied in Husebye et al. (1990) (see Section 2.4 for further description). The aim of their analysis is to estimate the mean length of the Migratory Motor Complex (MMC) period (i.e., the mean interoccurrence time). This data set was also analyzed in Aalen and Husebye (1991) using a variance component model and an intensity-based formulation with a gamma frailty component using a parametric Weibull model. Then Pe˜na et al. (2001) analyze this data using the estimators described in section 1.3. Although Aalen and Husebye (1991) stated that “the consecutive MMC periods for each individual appear (to be) approximate renewal process” we need to verify this assumption. To do so, Pe˜na et al. (2001) suggested that since formal statistical methods for checking this i.i.d. assumption are not yet available, a graphical method may be employed to assess the viability of the i.i.d. model by comparing the agreement among the PSH, WC, and FRMLE estimators in Pe˜na et al. (2001). The resulting estimates of the interoccurrence time survivor function are presented in Figure 3.2. A close agreement among these three estimates is evident. Thus, this agrement provides support for Aalen and Husebye’s assumption that the independence assumption is true. After showing that we can assume that the data follow an IID model, we then estimate survival function using either PSH, WC or FRMLE estimators. We can also estimate pointwise confidence intervals of the survival function which can be seen in the Figure 3.3. 3.4 Examples 57 50 100 150 200 250 0.0 0.2 0.4 0.6 0.8 1.0 Time Survivor Probability Estimates PSH (i.i.d.) WC FRMLE (α = 10.17) Figure 3.2: Survival function of interoccurrence times estimated using PSH, WC and FRMLE estimators for the MMC data set. This example corresponds to the IID model. 50 100 150 200 250 0.0 0.2 0.4 0.6 0.8 1.0 Time Survivor Probability Estimates Figure 3.3: Survival function of interoccurrence times estimated using PSH estimators and their pointwise 95% confidence interval for MMC data set using log-log transformation. 58 Chapter 3 Confidence Intervals for Median Survival 0 500 1000 1500 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Time Probability of rehospitalisation PSH (i.i.d.) WC FRMLE (α = 1.04) Figure 3.4: Probability distribution function of interoccurrence times estimated using PSH, WC and FRMLE estimators for the hospital readmission in colorectal cancer data set. This example corresponds to a frailty model with frailties following a gamma distribution. 3.4.2 Colorectal cancer rehospitalizations The next data set concerns rehospitalization of patients diagnosed with colorectal cancer analyzed in Gonz´alez et al. (2005). A description of this data set can be found in Section 2.1. As in the previous example we need to determine if the i.i.d. model is viable. The resulting estimates of the readmission time distribution are presented in Figure 3.4. We have displayed the estimates of the distribution function instead of the survival function because in this study, the investigator is interested in analyzing the probability of readmission instead of the probability of not visiting the hospital. A considerable difference between these three estimates is evident. The difference is clear between PSH and both WC and FRMLE estimators. Thus, basing on Pe˜na, Strawderman and Hollander’s argument, we may conclude that the i.i.d. model is not appropriate for this readmission data set. In making practical conclusions, it behooves therefore to use the inference obtained from the gamma frailty model. In addition, the estimate of the frailty parameter, α, is near 0 3.4 Examples 59 indicating a high correlation in the interoccurrence times. Thus, following these results, it was decided to analyze data using a Cox gamma shared frailty model (see Gonz´alez et al., 2005). Median survival times and probabilities to be readmitted at one and three years using FRMLE model can be found in Table 3.7. We realize that the one-year probability of rehospitalisation was 0.26 in Dukes stages A-B, 0.38 for patients with stage C, and 0.64 for those with stage D. This indicates that the probability of being readmitted strongly depends on advanced tumor stages. The three-year probability shows a similar difference. Males, chemotherapy use, high educational level, mortality, and high co-morbidity were also associated with smaller times between readmissions. After that, in order to verify if these observed differences in Dukes stage are statistically significant we compute their confidence intervals. Table 3.8 shows the confidence interval obtained using asymptotic variance and bootstrap procedures described in Section 3.1. As a general comment we can say that, and as we expected, that median survival time differs depending on the model. Under i.i.d assumption median survival is underestimated. Regarding Dukes stage, we can conclude that there are differences between median survival time for patients with stage D and both patients with stages A-B and C. These differences are statistically significant since their confidence intervals does not overlap. This conclusion is the same assuming both independent and correlated model. We also observe that the αestimates are very small in all cases (Table 3.8). So, it seems reasonable to use a frailty model for making these comparisons. Regarding differences in confidence intervals and their width, we should mention that these results agree with those observed in the simulation study. We first consider the confidence intervals computed using asymptotic variances. Table 3.8 shows that the narrowest confidence intervals are those obtained using the log-log transformation (joint with arcsinus in some occasions) as simulation studies showed (Tables 3.5 and 3.6). On the other hand, when bootstrap method is used, confidence intervals are close to those obtained using log-log transformation although they are slightly wider, as the simulations also indicated. 66 Chapter 3 Confidence Intervals for Median Survival Needs to Determine a Seed Value for Alpha Seed Alpha: 12.55364 Alpha estimate= 1.460141 Needs to Determine a Seed Value for Alpha Seed Alpha: 9.342046 Alpha estimate= 2.193977 # Median survival estimates > print(fit.FRMLE,digits=c(4,2)) Survival for recurrent event data. Group= as.factor(dukes) n events mean se(mean) median recurrences: min max median 1 180 144 1350.0 67.50 2175 0 6 0 2 148 183 841.5 46.82 1073 0 16 1 3 75 131 363.2 45.96 199 0 22 1 # Probability at one and three years > for (i in 1:length(fit.FRMLE)) + { + cat(i,"\n") + cat("at one year: ",1-fit.FRMLE[[i]]$surv[sum(fit.FRMLE[[i]]$time<365)],"\n") + cat("at three years: ",1-fit.FRMLE[[i]]$surv[sum(fit.FRMLE[[i]]$time<1095)],"\n") + } 1 at one year: 0.2604853 at three years: 0.4072520 2 at one year: 0.3814841 at three years: 0.5045216 3 at one year:0.6418163 at three years: 0.8867255 Asymptotic confidence intervals, for WC estimator, described in Section 3.1 can be obtained in Ras follows: > fit.WC<-survfitr(Survr(id,time,event)~as.factor(dukes),data=readmission,t="w") > cat("dukes=1 \n") dukes=1 3.5 R instructions for survrec package 67 > print(Brook.Crowley(fit.WC$"1",0.5)) $percentile [1] 1736 $ci95.asymptotic $ci95.asymptotic$bandwith [1] 1148 $ci95.asymptotic$ci95 [1] 1446.15 2025.85 $ci95.id [1] 655 1736 # formula 3.8 $ci95.log.log [1] 1157 1736 # formula 3.9 $ci95.arcsin [1] 1157 1736 # formula 3.1 > cat("dukes=2 \n") dukes=2 > print(Brook.Crowley(fit.WC$"2",0.5)) $percentile [1] 1028 $ci95.asymptotic $ci95.asymptotic$bandwith [1] 735 $ci95.asymptotic$ci95 [1] 589.27 1466.73 $ci95.id [1] 276 1483 # formula 3.8 $ci95.log.log [1] 489 1291 # formula 3.9 $ci95.arcsin [1] 462 1291 # formula 3.1 > cat("dukes=3 \n") dukes=3 > print(Brook.Crowley(fit.WC$"3",0.5)) $percentile [1] 199 $ci95.asymptotic $ci95.asymptotic$bandwith [1] 132 $ci95.asymptotic$ci95 [1] 116.87 281.13 $ci95.id [1] 79 474 # formula 3.8 $ci95.log.log [1] 158 297 # formula 3.9 $ci95.arcsin [1] 158 297 # formula 3.1 where Brook.Crowley function has been created because it is not included in the survrec package (see Appendix E) On the other hand, the same confidence intervals using bootstrap procedures described in Section 3.2 can be fitted using survdiffr function. Let us compute bootstrap Plan II, that is, nonparametric bootstrap estimating ˆ Fusing WC estimator, and simulating Gfrom its empirical distribution. 68 Chapter 3 Confidence Intervals for Median Survival > fit<-survdiffr(Survr(id,time,event)~as.factor(dukes),data=readmission, + q=0.5,boot.F="WC",boot.G="empirical",B=999) > print(fit) $"1" CASE RESAMPLING BOOTSTRAP FOR CENSORED DATA Call: survdiffr(formula = Survr(id, time, event) ~ as.factor(dukes), data = readmission, q = 0.5, B = 999, boot.F = "WC", boot.G = "empirical") Bootstrap Statistics : original bias std. error t1* 1736 -459.972 810.8057 $"2" CASE RESAMPLING BOOTSTRAP FOR CENSORED DATA Call: survdiffr(formula = Survr(id, time, event) ~ as.factor(dukes), data = readmission, q = 0.5, B = 999, boot.F = "WC", boot.G = "empirical") Bootstrap Statistics : original bias std. error t1* 1028 -54.57057 251.8067 $"3" CASE RESAMPLING BOOTSTRAP FOR CENSORED DATA Call: survdiffr(formula = Survr(id, time, event) ~ as.factor(dukes), data = readmission, q = 0.5, B = 999, boot.F = "WC", boot.G = "empirical") Bootstrap Statistics : original bias std. error t1* 199 14.52452 46.15239 We notice that survdiffr function returns an object of class "boot". Thus, we can use boot package for summarizing the object fit using generic print function. Then, bootstrap confidence intervals showed in Table 3.8 can be fitted using boot.ci. This function is also included in boot package. > # Dukes stage A-B > boot.ci(fit$"1") BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 999 bootstrap replicates 3.5 R instructions for survrec package 69 CALL : boot.ci(boot.out = fit$"1") Intervals : Level Normal Basic 95% (1148, 2347 ) (1297, 2315 ) Level Percentile BCa 95% (1157, 2175 ) (1188, 2175 ) Calculations and Intervals on Original Scale > # Dukes stage C > boot.ci(fit$"2") BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 999 bootstrap replicates CALL : boot.ci(boot.out = fit$"2") Intervals : Level Normal Basic 95% ( 608, 1566 ) ( 731, 1567 ) Level Percentile BCa 95% ( 489, 1325 ) ( 453, 1325 ) Calculations and Intervals on Original Scale > # Dukes stage D > boot.ci(fit$"3") BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS Based on 999 bootstrap replicates CALL : boot.ci(boot.out = fit$"3") Intervals : Level Normal Basic 95% ( 95.3, 275.3 ) ( 48.0, 240.0 ) Level Percentile BCa 95% (158, 350 ) (113, 247 ) Calculations and Intervals on Original Scale Some BCa intervals may be unstable Chapter 4 Inference for the General Class of Models This chapter gives procedures for estimating the parameters of the general class of models for recurrent events proposed by Pe˜na and Hollander (2004). An approach based on semiparametric inference (Pe˜na, Slate, and Gonz´alez, 2003) has been developed. When model without frailties is estimated, a generalization of the partial likelihood for the Cox model is obtained. On the other hand, when frailties are included in the model, an EM algorithm is developed. As we have mentioned in Section 1.6, this work was not only carried out by myself, so we have included it in the Appendix B. We encourage the reading of this appendix before going through this chapter because the present chapter is based on the notation and the results we get in it. Herein, we present another methodology based on a penalized likelihood approach. Two different strategies are adopted. One of them was developed in the shared frailty model context by Therneau et al. (2003). Their idea is based on penalizing the partial likelihood where the penalization bears on a regression coefficient (see Section 1.4.2 for further details). The second penalized approach, also applied in the shared frailty model, was proposed by Rondeau et al. (2003). Their method of estimation is based on the penalized full likelihood, and it gives a nonparametric estimation of the baseline hazard function using a continuous estimator. The solution is then approximated using splines. The motivations for the present chapter are mainly due to general drawbacks which appears when EM algorithm is used. In particular, direct estimates of the variance of parameters are 71 72 Chapter 4 Inference for the General Class of Models not provided. Thus, to solve this problem, we first propose to estimate parameters involved in Pe˜na and Hollander’s model using the ideas described in Therneau et al. (2003). However, this approach still continue having problems such as the convergence can be slow and the variance of frailty variance cannot be estimated directly. To improve this approach, an alternative method of penalization proposed by Rondeau et al. (2003) is adopted. The main advantage of this method is that we can easily obtain smooth estimates of the hazard function and an estimation of the variance of frailty variance parameter. In Section 4.1 we begin by giving how to fit Pe˜na and Hollander model using EM algorithm. We also compare the resulting estimates using this approach with those obtained using existing methods such as AG, WLW and PWP models described in Section 1.4.3. Section 4.2 deals with both methods of estimation based on penalization above outlined. In Section 4.3, based on Nielsen et al.’s Nielsen et al. work, we give some procedures to test whether the frailty is necessary. Some computational issues are discussed in Section 4.4. In Section 4.5 the estimation procedures are illustrated using both readmission and bladder cancer data sets. This last example is also used to compare the results obtained using Pe˜na and Hollander’s model to those obtained using three methods of analyzing recurrent event data mentioned in Chapter 1. Bladder cancer data set is also used to illustrate which is the impact of miss-specifying effective age process. Finally, Section 4.6 shows how to fit the general class of models using gcmrec and frailtypack packages. 4.1 Semiparametric inference: EM algorithm Before going through the penalized likelihood inference we focus on presenting how to fit Pe˜na and Hollander model using EM algorithm described in Appendix B. To do this, we consider the data set belonging to patients with colorectal cancer described in Section 2.1. In this analysis we consider only the variables tumor stage (Dukes classification: A-B, C or D) and gender. This example will also be useful to describe how to fit and interpret the results obtained using the general class of models as well as compare them to the typical Cox models for recurrent event data described in Section 1.4.3. This data set was first analyzed in Gonzalez et al. (2005) using a shared gamma frailty model since interoccurrence times were correlated. Thus, instead of using a marginal models such as those described in 1.4.3 we have fitted both AG and PWP models including a frailty term. Table 4.1 Semiparametric inference: EM algorithm 73 4.1 shows the resulting estimates using these approaches together with those obtained using the general class of models. First of all, we observe that there is a significant random effect since ξ (frailty precision, see Appendix B) is quite small in all cases. When the effective age corresponds to perfect repair, the resulting estimates from Pe˜na and Hollander model are close to those obtained with the PWP conditional method. On the other hand, when a minimal repair effective age formulation is used, the results are close to those obtained using the AG model although the hazard risk are little different. The differences in hazard estimates are probably due to the impact of accumulating event readmissions. This effect is incorporated in the Pe˜na and Hollander model via the ρfunction. Using, ρ(k;α) = αkwe model different scenarios. As an example, if αis less than unity, the increasing number of rehospitalization has a beneficial effect. In our case, the probability of being rehospitalized would decrease with the number of hospitalization. Looking at our results (Table 4.1) we observe that αis greater than unity, indicating that each hospitalization increases the risk of further hospitalization. Using other words, αparameter greater than one suggests that there is different risk of being hospitalized depending on the number of rehospitalizations. We may further analyze this fact looking at the cumulative hazard functions for the time since last event estimated using the PWP stratified model. The resulting plot, shown in Figure 4.1, shows that the probability of being readmitted increases as the number of hospitalization increases (red lines). On the other hand, the cumulative hazard estimates using Pe˜na and Hollander model is common for all events, and approximates to the average of all marginal cumulative hazards. For this model, the alpha parameter allows changing the baseline hazard and it play the same role as the stratification in the PWP model. Regarding risk estimate, as PWP model is an stratified model, we can fit separate coefficients for both sex and dukes variables to each stratum. Table 4.2 contains these coefficients where the rehospitalizations greater than fourth are combined. We can observe that the differences between males and females are only statistically significant in the probability of being readmitted for the second rehospitalization. This fact is difficult to be explained as the differences observed between males and females using a simple fit model (Table 4.1) as Gonzalez et al. (2005) pointed out. On the other hand, patients with advanced tumoral stage (dukes D) have and approximately 4.5-fold first recurrence rate as compared to those with early tumoral stage (dukes A-B). This difference decrease when posterior rehospitalizations are analyzed to be approximately 2-fold risk. Patients 74 Chapter 4 Inference for the General Class of Models Pe˜na and Hollander model Shared Gamma Frailty model perfect repairsaminimal repairsbPWP model AG model Covariate HR (CI95%) HR (CI95%) HR (CI95%) HR (CI95%) Gender Female 1 1 1 1 Male 1.56 (1.20-2.03) 1.69 (1.19-2.40) 1.47 (1.15-1.88) 1.85 (1.34-2.54) Dukes stage A-B 1 1 1 1 C 1.49 (1.08-2.05) 1.57 (1.04-2.36) 1.46 (1.12-1.91) 1.63 (1.15-2.31) D 3.06 (2.04-4.58) 4.05 (2.25-7.29) 3.45 (2.47-4.81) 5.03 (3.35-7.56) Frailty ξ2.50 1.03 3.09 0.73 log N(s−)α1.08 (0.97-1.19) 1.11 (0.61-1.60) aEffective Age is backward recurrence time (E(s) = s−SN†(s−)) bEffective Age is calendar time (E(s) = s). Table 4.1: Hazard ratios and 95% confidence intervals for the probability of rehospitalization to the readmission data set. Estimates using Andersen-Gill (AG) and Prentice, Williams and Peterson (PWP) methods with a frailty term, together the estimates obtained using Pe˜na and Hollander model using two different effective age formulations. with dukes C only differ from patients with early tumoral stages (dukes A-B) in the probability of having the first hospitalization (HR:1.66, CI95%: 1.19-2.75). 4.2 Penalized Likelihood Inference Herein, we propose a different method of parameter estimation that solves problems which appear when EM algorithm is used. In Section 4.2.1, we begin by showing how to adapt Therneau et al.’s 2003 approach in our case while in Section 4.2.2 we develop another different penalized approach that is able to give an estimation of the variance of frailty variance parameter. This gives us a possibility making a formal test to see whether data are correlated or not. Let us assume that Xiis time-independent, an increasing number of event occurrences is of form ρ(N† i(s−); α) = αN† i(s−)(that is, αk, where kis the number of occurrence), and ψ(x) = exp(x). Assuming this situation, the Pe˜na and Hollander model (1.20) can be also written as follows: λi(s|Zi,Xi) = Ziλ0[Ei(s)] αN† i(s−)exp[β0Xi],(4.1) As in Pe˜na et al. (2003), we also assume that the frailties Z1, . . . , Znare i.i.d. from a gamma distribution. However, in our case we consider the parameterization ν= 1/ξ. Therefore, we 4.2 Penalized Likelihood Inference 75 Covariate HR (CI95%) p value Male 1st event 1.26 (0.93-1.70) 0.130 2nd event 1.80 (1.16-2.79) 0.009 3rd event 0.92 (0.49-1.72) 0.790 ≥4 1.60 (0.78-3.28) 0.200 Dukes C 1st event 1.66 (1.19-2.31) 0.003 2nd event 1.36 (0.83-2.24) 0.220 3rd event 0.82 (0.39-1.72) 0.600 ≥4 1.20 (0.53-2.75) 0.660 Dukes D 1st event 4.54 (3.12-6.59) <0.001 2nd event 2.15 (1.22-3.78) 0.008 3rd event 2.16 (1.02-4.56) 0.043 ≥4 2.03 (0.82-5.03) 0.130 Table 4.2: Hazard ratios and 95% confidence intervals per event for the probability of rehospitalization to the readmission data set. Estimates using stratified Prentice, Williams and Peterson (PWP) model. 0 500 1000 1500 0.0 0.5 1.0 1.5 Days since last event Cumulative Hazard PWP model event 1 event 2 event 3 event >=4 Pena and Hollander Figure 4.1: Baseline cumulative hazards for each event number using Prentice, Williams and Peterson (PWP) conditional method and Pe˜na and Hollander model to readmission data set. 82 Chapter 4 Inference for the General Class of Models 4.4 Computational Issues The estimation procedures for the general class of models proposed by Pe˜na and Hollander (2004) have been implemented in an Rpackage which is available at the CRAN project (Ihaka and Gentleman, 1996, R Development Core Team, 2005, URL:http://www.R-project.org). The function gcmrec performs estimation for model (1.20) with or without frailties, with ρ(j;α) = 1 or ρ(j;α) = αj, and with ψ(·) = exp(·). These procedures are implemented using the combination of R(flexible, high-level statistical language) with Fortran (high execution speed of iterative procedures). Thus, we have implemented a dynamic link library (dll) in Fortran 77 that is called by gcmrec Rfunction. For efficiency, the numerically intensive steps of the algorithm are coded in Fortran 77, loaded as a dynamic linked library into R, and invoked from the top-level Rroutine (also named gcmrec). A major problem in the parameter estimation is the maximization procedures of the likelihood. We comment briefly on some aspects of methods used in the programming. Estimation of αand βrequires maximizing the logarithm of the profile likelihood in the model without frailty, and, additionally in the model with frailty, obtaining the conditional expectation of {Zi}. Newton-Raphson (N-R) algorithm can be used to obtain solutions as well as to obtain the inverse of the approximate Hessian. The algorithm continues until convergence, or until a pre-set maximum number of iterations is reached. This algorithm may diverge if the Hessian is not positive definite. To circumvent this, we can modify the Hessian by adding a large enough constant to its diagonal, making the matrix positive definite but perturbing the Hessian as little as possible. This modification requires computing the Choleski factorization in some steps of the N-R procedure. Both the N-R and modified Hessian procedures are implemented in gcmrec function. When we perform estimation with frailty models by penalizing partial likelihood, memory and time-consuming considerations may become an important issue. For instance, let us assume that we have 403 patients such as in the hospital readmissions data set and that we fit the model with frailties. Let us also assume that we are including 6 other variables. Then, the full information matrix has 4092= 167281 elements. One possible solution to this problem is to use quasi-Newton methods that can be used when the Hessian matrix is difficult or time-consuming to evaluate. In this case, instead of obtaining an estimate of the Hessian matrix at a single point, these methods gradually build up an approximate Hessian matrix by using gradient information from some or 4.4 Computational Issues 83 all of the previous iterations. BFGS algorithm is one of the widely used quasi-Newton methods. Another possibility is to use a modification of this algorithm, which is known as L-BFGS-B that is a limited-memory quasi-Newton code for large-scale bound-constrained or unconstrained optimization (see Byrd and Nocedal, 1995 for further details). The main disadvantage of these algorithms is that as the inverse Hessian is estimated at each step, we cannot use it for estimating the variance of the parameters. Thus, in our programs, we adopted another approach that was proposed by Therneau et al. (2003) in the case of shared frailty model. Their idea was to partition the inverse of the Hessian according to the rows of Xand M(that makes reference to covariates and frailties, respectively) as follows: I= IMM IMX IXM IXX  . Then, they propose to use a sparse computation option, where only the diagonal of IMM is retained. They indicate that this method has not a large impact on the estimation procedure because neither the score vector nor the likelihood are changed. Thus, the solution is identical to the one obtained in the non-sparse case, but the speed of the algorithm increased dramatically (see Therneau et al., 2003 for further details). Finally, when penalized full likelihood approach is used, both the first derivative (the score) and the second derivative (the Hessian) of the log-likelihood themselves do not have a simple analytical form, so another different method from N-R must be used. To solve this problem, we have chosen the Marquardt algorithm which computes numerically these quantities using finite differences. This algorithm was first published by Levenberg (1944). Then, it was rediscovered by Marquardt (1963) who applied it to statistical problems. The Levenberg-Marquardt algorithm (LMA) provides a numerical solution to the problem of maximizing a sum of square of several, generally nonlinear functions that depend on a common set of parameters. This method is applied in our case by maximizing the sum of squares of the partial derivatives to find their solutions. The LMA is a combination between the N-R algorithm and the method of gradient descent and is more robust than the N-R, which means that in many cases it finds a solution even if it starts very far off the final minimum. In our problem, to be sure of having a positive function at all stages of the algorithm, we restrict all the spline coefficients ηjto be positive for all j. We imposed a constraint of positivity for the parameter ν, so we did not consider a negative dependence in the model, which obviously 84 Chapter 4 Inference for the General Class of Models do not have a frailty interpretation (although negative values make sense in the intensity as we have illustrated in Section 4.3). In our approach we use a modified Marquardt algorithm like the Newton-Raphson procedure. Here, we outline the main points of this algorithm also used in Rondeau et al. (2003). Let θbe the parameters to be estimated (in our case (η, α, β, ν)). If necessary the diagonal of the Hessian at iteration k,H(k), is inflated to obtain a positive definite matrix H∗(k)=H∗(k) ij with H∗(k) ii =H(k) ii +%[(1 −δ)|Hk ii|+δtrace(H)] where %and δare parameters with initial values set equal to 0.01. They are reduced when H∗is positive definite and increased if not, and H∗(k) ij =H(k) ij . The estimates θ(k)are then updated to θ(k+1) using both the current modified hessian, H∗(k), and the current gradient of the parameters according to the formula: θ(k+1) =θ(k)−φ∇(θ(k)) H∗(k) where ∇denotes the gradient and if necessary, φis modified to ensure that the log likelihood is improved at each iteration. One-dimensional maximization In the model with frailties, we further need to maximize the marginal likelihood with respect to only one parameter, ν. The gcmrec package provides for two options here: the Newton-Raphson method, and Brent’s algorithm (Brent, 1973) which has a faster linear rate of convergence than golden section search. This method is as a one-dimensional maximization without derivatives. First we bracket the maximizing value, and then we obtain it using Brent’s method in one dimension (see Brent, 1973 for further details). In both cases, optimization is performed using the reparameterization to log(ν). Another possibility is to use ν1/2as Therneau and Grambsch (2000) used in the context of shared frailty model. Test of Rfunctions Therneau and Grambsch (2000, Appendix E) gave a set of test data with known answers. These data sets were mainly created to illustrate the computations of statistical algorithms and to ensure 4.5 Hospital Readmission and Bladder Cancer Data Sets Revisited 85 the accuracy and quality of software programs (applied in the survival analysis settings). As the Cox model and shared frailty model are particular cases of the Pe˜na and Hollander model, we have used these data sets to test gcmrec function. Using our function, we reproduce the exact results that are obtained using survival library (that used in Therneau and Grambsch, 2000). For the Cox model we set up Z= 1, ρN† i=1 and Ei(s) = s−SiN† i(s−)in the gcmrec function and we obtain the same results as using coxph function. Similarly, the shared frailty model (using coxph and frailty functions in survival package) showed the same results as using gcmrec function when ρN† i=1 and Ei(s) = s−SiN† i(s−). 4.5 Hospital Readmission and Bladder Cancer Data Sets Revisited In this section we apply the estimation procedures developed in preceding sections to two real data sets previously analyzed in Chapter 3 using a more complex model. 4.5.1 Hospital Readmission Study In this example we consider only the variables tumor stage (Dukes classification: A-B, C or D) and gender corresponding to the data about rehospitalization in patients with colorectal cancer described in Section 2.1. Since in this example we do not have information about the effective age, we assume the backward recurrence time, E(s) = s−SN†(s−). We fitted the general model (with frailties), taking s∗= 2060, the maximum follow-up time. After 31 iterations in the EM algorithm (see Section 4.6), the estimate of the frailty precision is ξ= 1/ν is quite small (ˆ ξ= 2.50), so we may conclude that the frailty component of the model is important for these data. We also saw this fact in Section 3.4.2 by using a graphical method. However, we cannot conclude its statistical signification since the EM algorithm does not provide any estimation of frailty variance. On the other hand, we may also say that among these covariates, the advanced tumor stages (C or D) and males are associated with an elevated risk of rehospitalization. Furthermore, since the estimate of αis larger than unity, there is an indication that each hospitalization increases the risk of further hospitalization, as one could expect. We may use another method of estimation such as penalized likelihood from the two different points of view we have outlined in Section 4.2. Table 4.3 shows the hazard ratios and its confidence 86 Chapter 4 Inference for the General Class of Models Penalized approach EM approach Partial Likelihood Full Likelihood Jackknife Covariate HR (95%CI) HR (95%CI) HR (95%CI95) Gender Female 1 1 1 Male 1.56 (1.17-2.09) 1.63 (1.25-2.13) 1.56 (1.20-2.03) Dukes stage A-B 1 1 1 C 1.49 (1.11-2.00) 1.51 (1.14-2.01) 1.49 (1.08-2.05) D 3.06 (2.14-4.39) 3.41 (2.39-4.86) 3.06 (2.04-4.58) log N(s-) α1.08 (0.97-1.19) 1.14 (1.03-1.24) 1.08 (0.81-1.35) Frailty ν0.40 0.58 0.40 (SE ν) (NA) (0.15) (NA) ξ2.50 1.72 2.50 κ3.36 ×1011 Table 4.3: Hazard ratios (HR) and 95% confidence intervals for the probability of rehospitalization for the colorectal data set. Estimates using both penalized and EM approach assuming effective age gap time formulation, E(s) = s−SN†(s−). Standard errors are computed using H−1IH−1for penalized approach and using Jackknife for the EM approach. 4.5 Hospital Readmission and Bladder Cancer Data Sets Revisited 87 intervals using the EM approach and both methods of penalization. First of all we deal with the importance of frailty. As we have previously mention, we have used a graphical method to check independence assumption. Now, using the penalized full likelihood approach, standard error of the frailty variance can be estimated. So, we can then verify independence assumption using one-side Wald statistic (see Nielsen et al., 1992 or Self and Liang, 1987 and also Section 4.3). In the model ˆν/SE(ˆν) = 0.58/0.15 = 3.87 and 1-pnorm(3.87)=5.441768e-05. Another possibility to test the null hypothesis H0:ν= 0 is to use a likelihood ratio test as Therneau and Grambsch (2000) or Nielsen et al. (1992) indicated. As we have mentioned in Section 4.3, this likelihood ratio test may be performed as twice the difference between the log-partial-likelihood with the frailty terms integrated out, and the log-likelihood of a model without frailties. In our case, the test for significance of the frailty is −2751.2 vs. −2719.7, which gives a chi-square statistic of 31.5 on one degree of freedom for a p-value of 2.11e-15. So, one can conclude that there is heterogeneity among the interoccurrence times. We notice that that EM algorithm and penalized partial likelihood approach should be identical as Therneau and Grambsch (2000) proved in the case of shared frailty model. The prove of this results for the Pe˜na and Hollander model is beyond the scope. We also realize that there is some little differences between the risks estimates using EM algorithm and penalized full likelihood, probably due to differences in the maximization procedures and the selection of number of knots and bandwidth. However, in all case we can conclude that males and patients with advance tumoral stage have more probability to be readmitted and that these risks are statistically significant. On the other hand, we also can see as the width of confidence intervals for hazard ratios are very different depending on the method of estimation used. Having observed these differences, it is of interest to compare the estimates of variance of α and βparameters using the three different approaches. To do so, we have carried out a little simulation study following Pe˜na et al.’s 2003 procedure (this simulation study is also described in Section 5.3). Table 4.4 summarizes the standard deviation estimates using EM algorithm and both penalized likelihood approaches: partial (PPL) and full (PFL). From this table we note that EM algorithm and PPL approach overestimate the empirical standard deviation of parameters associated with the covariates, β, while the variance of αestimates is quite well estimated. On the other hand, PPL approach clearly underestimates the empirical standard deviations, specially for the αparameter. Thus, after our results, we can say that the conclusions obtained after analyzing 88 Chapter 4 Inference for the General Class of Models γ ξ ˆσˆαˆσJACK ˆαˆσPPL ˆαˆσPFL ˆαˆσˆ β1ˆσJACK ˆ β1ˆσPPL ˆ β1ˆσPFL ˆ β1ˆσˆ β2ˆσJACK ˆ β2ˆσPPL ˆ β2ˆσPFL ˆ β2 0.9 2 0.007 0.008 0.007 0.003 0.387 0.562 0.483 0.331 0.228 0.312 0.293 0.196 2 2 0.006 0.007 0.007 0.002 0.352 0.498 0.475 0.326 0.180 0.271 0.261 0.166 0.9 6 0.007 0.007 0.006 0.002 0.262 0.403 0.382 0.244 0.162 0.241 0.238 0.149 2 6 0.006 0.007 0.008 0.002 0.221 0.331 0.341 0.204 0.129 0.196 0.185 0.101 0.9 ∞0.007 0.007 0.008 0.002 0.163 0.221 0.214 0.159 0.108 0.183 0.168 0.081 2∞0.005 0.006 0.006 0.001 0.128 0.202 0.196 0.131 0.074 0.158 0.146 0.056 Table 4.4: Summary of empirical standard deviations of the estimators of α(ˆσˆα), and β(ˆσˆ β1and ˆσˆ β2) using EM algorithm and penalized likelihood approaches (JACK: jackniffe estimates, PPL: penalized partial likelihood, PFL: penalized full likelihood). This corresponds to the simulation study performed in Pe˜na et al. (2003) for the case where the αparameter is 1.05 and the sample size n= 30. The true value of βis (1,−1), and 1000 replications were run for each parameter combination. The other parameter are: γthe Weibull shape and ξthe frailty precision. readmission data set (Table 4.1) are in general correct, except for the statistically significance of αparameter stem from PFL approach. The last analysis performed in this data set was to compare the baseline survivor function estimates using EM algorithm and PFL approach. Figure 4.2 shows a graphical comparison between the estimation of baseline survivor function using both approaches. We can see that the curve provided by PFL method is completely smooth, while EM algorithm gives a function which jumps at each observed readmission time. 4.5.2 Bladder Cancer Study We analyze the covariates: X1, the treatment indicator (1 = placebo; 2 = thiotepa); X2, the size (in cm) of the largest initial tumor; and X3, the number of initial tumors. First, we fit the Pe˜na and Hollander model using the gap time formulation, E(s) = s−SN†(s−), as effective age. With s∗= 64, the maximum observation period, the general model without frailties estimates: ˆα= 0.9826 and (ˆ β1,ˆ β2,ˆ β3) = (−0.3188,−0.0154,0.1353). Same results are also obtained using model with frailties since in that case ˆ ξ= 5432999 (ˆν≈0). Thus, using the approximate inverse of the partial likelihood information matrix from fitting the model without frailties, the associated estimated standard errors are .0736 for ˆαand (0.2051,0.0695,0.0511) for ˆ β(Table 4.5). The effective age for these data is not known, so we also fitted the general model with frailties assuming a calendar time for effective age, E(s) = s. For this case, the estimates are ˆα=.789, (ˆ β1,ˆ β2,ˆ β3) = (−.5743,−.0315, .2220), and ˆ ξ=.974 (ˆν= 1.03). Here we could say that the frailty 4.5 Hospital Readmission and Bladder Cancer Data Sets Revisited 89 0 500 1000 1500 2000 0.0 0.2 0.4 0.6 0.8 1.0 Time Baseline survivor function EM algorithm Penalized full likelihood Figure 4.2: Baseline survivor function estimated using both EM algorithm and penalized full likelihood approach for time until next rehospitalization corresponding to readmission data set. parameter is important. In order to compare both models, the estimates of the survivor functions for the two effective age specifications are presented in Figure 4.3. The lower curves (blue lines), corresponding to the placebo group, are obtained by setting X1= 1 in the expression given by {ˆ ¯ F0(t)}exp{ˆ β1X1+ˆ β2¯ X2+ˆ β3¯ X3}, while the upper curves (dark green lines) are for the thiotepa group obtained by setting X1= 2. The observed means were ¯ X2= 2.01 and ¯ X3= 2.11. The solid curves are for the backward recurrence time effective age, while the dashed curves are for E(s) = s. We observe that thiotepa group shows a higher survival rate than the placebo group, although the statistical significance of this difference depends on which effective age process was used. Then, we compare these results with those obtained using the three existing methods of analysis described in Therneau and Hamilton (1997) and Therneau and Grambsch (2000). Table 4.5 summarizes the estimates from, AG, WLW, and PWP methods, together with the estimates obtained from the Pe˜na and Hollander model using two specifications of the effective age process, E(s) = s−SN†(s−)and E(s) = s. The authors analyzed the data set by using models described in Section 1.4.3 (AG,PWP and WLW). These models are marginal models, while the general class 90 Chapter 4 Inference for the General Class of Models 0 10 20 30 0.0 0.2 0.4 0.6 0.8 1.0 Time Survivor Function Perfect repair (gap time) Minimal repair (calendar time) Figure 4.3: Estimates of the survivor function for bladder data set when the Pe˜na and Hollander model is fitted. The blue curve corresponds to the placebo group, while the dark green curve is for the thiotepa group, both evaluated at the mean values of size of initial tumor and mean number of initial tumors. The solid curves show effective age E(s) = s−SN†(s−)(perfect repair), while the dashed curves are when E(s) = s(minimal repair). 4.5 Hospital Readmission and Bladder Cancer Data Sets Revisited 91 of models is a frailty model. As we have mentioned in Section 1.4.2, the estimates from a frailty model have a subject-specific interpretation, while marginal models have a population-average interpretation. Hence, the estimates are not directly comparable to those presented in Therneau and Hamilton (1997) or in Therneau and Grambsch (2000). To make the results comparable we fitted both AG and PWP models including a frailty term (Table 4.5). First of all, we highlight the role that the effective age process plays in data analysis. We also note that the results obtained using some of existing models are similar to those obtained using the general model. As an example, when ‘perfect repairs’ is assumed, the results obtained are close to those obtained using PWP model. In both cases the frailty term was not important. On the other hand, the results obtained using ‘minimal repairs’ are close to those obtained from both AG with frailties and WLW method. We see that in this case it is necessary to include a frailty term to model the association among the inter-event times for each patient. This similitude is what one would expect since the time scale acting in the hazard for the PWP model is gap time, E(s) = s−SN†(s−), while in the case of WLW method is time scale formulation, E(s) = s. We realize that the parameter estimates using the AG model without covariates lie between WLW and Pe˜na and Hollander model (estimates from Therneau and Grambsch, 2000) indicating the need of using a frailty model. The observed differences between models indicate the importance of the effective age and the need to monitor this information. This example along with the results we are showing in the next section, are the basis to further study how to incorporate the notion of effective age in biomedical setting. This was our goal in Chapter 5 where we illustrate how to incorporate information about intervention after relapses in patients with cancer. 4.5.3 Miss-specification of effective age Because of the importance of the effective age process as demonstrated by previous application to the bladder cancer data, we examined further through a simulation study the impact of missspecifying the effective age process. We again consider the simulation study carried out in Pe˜na et al. (2003). We examine the impact of two types of effective age process miss-specification: that the interventions following event occurrences are all minimal repair, or that they are all perfect repair. Tables 4.6 and 4.7 show the results for a given simulation (when α=.9). The results indicate an interesting interplay between the nature of the baseline survivor func- 98 Chapter 5 Dynamic Cancer Model for Tumor Relapses related to the time until progression is important because most patients die from causes related to the disease (see Lister, 1991 or Romaguera et al., 1991). Thus, the estimation of the risk of recurrence would allow for better planning of follow-up after diagnosis or first treatment, and would permit clinicians to consider new therapeutic approaches depending on the patient’s risk of relapse. In Chapter 1 we have shown that there exists many survival models that handle recurrent event data. However, several prognostic studies in major cancer and epidemiologic journals estimate the risk of relapse only using information about the time until first occurrence (follicular lymphomas Lombardo et al. 2002, Lopez-Guillermo et al., 2000, acute leukemias Godder et al., 2004, colorectal cancer Schwandner et al., 2000, or breast cancer Fredriksson et al., 2002, among others). This approach ignores the information of subsequent relapses, hence statistical inference will tend to be inefficient. To avoid this problem, other cancer studies, such us mammary tumors for rats (see Gail et al., 1980) or patients with superficial bladder cancer (see Cheuvart, 1988 or Byar, 1980) rely on the Cox’s proportional hazards models and its variants which handle both intra-subject correlation and event dependence. As we have illustrated in Section 1.3.1, another aspect that may modify the event occurrence intensity arises from the interventions performed on the subject after each event occurrence. In cancer settings, patients with the disease are treated after observing a progression of the tumor. In the particular case of indolent lymphomas, as we are dealing with non-curable disease, the therapy aims to increase as much as possible the time until the next relapse. After giving some therapy, the patient is monitored and then we observe if cancer or some disease-related symptoms disappear. Thus, patients whose disease completely disappears will have less probability to relapse than those where little or no response is observed. Hence such interventions can be viewed as improving the patient. In the reliability literature this is referred to as adjusting the effective age of the system. This is the basis of work we developed in present chapter. We briefly outline the contents of this chapter. We begin by providing some additional notation for the general class of models in Section 5.1. This section also outlines how this class of models subsumes some existing models for recurrent event data which have been used in biomedical settings. Section 5.2 presents a description of the effective age process in cancer settings. A simulation study is performed to study the model behavior under effective age mis-specification. The simulation results and design are discussed in Section 5.3. Section 5.4 illustrates the use of 5.1 The Pe˜ na and Hollander Model Revisited 99 the model with a real data set from low grade lymphomas. Finally, Section 5.5 shows how to fit the dynamic cancer model using gcmrec package. 5.1 The Pe˜na and Hollander Model Revisited For a patient, as in the general case, in cancer settings we can observe a, possibly time-varying, q-dimensional vector of covariates such as gender, age, race, disease status, beta-2 microglobulin level, treatment regimen, etc. We denote this covariate process by {X(s) = (X1(s), X2(s), . . . , Xq(s))0: 0 ≤s≤τ}. In addition, in our case, after treatments or interventions are administered upon relapses, information about patient status may be obtained. Examples of interventions are chemotherapy, radiotherapy, and bone marrow transplant, among others. We denote this information by a vector ψ= (ψ1, ψ2, . . . , ψK)0, where ψjsignifies a certain type of response to the intervention after the jth relapse. This will be explained in more detail later. Consequently, if in the study there are npatients, over the period [0, τ∗] where τ∗≡maxi≤nτi, we will have the following data: D(τ∗)≡ {[(Xi(s) : 0 ≤s≤τ), ψi, Ki, τ, Si1, Si2, . . . , SiKi, τi−SiKi], i = 1,2, . . . , n}. We recall that the conditional intensity function given in 4.1 is λi(s|Z,Xi) = Ziλ0[Ei(s)] αN† i(s−)exp[βtXi(s)]. This model incorporates the effect of performed interventions through the effective age, Ei(s), which serves as the argument to the baseline hazard rate function Pe˜na and Hollander (2004). Here, we briefly show as some of existing models that have been used in recurrent events problems are particular cases of this general class of models. The simplest model regarding the effective age is when Z= 1, ρ[N† i(s−); α] = 1, and the effective age is always the same type and is either of the following possibilities: First, a patient can achieve a perfect response to the treatment. This means that the patient recovers perfectly and the status is the same as at the beginning of the study (if there is no time-dependent covariate). This model has been considered by Prentice 100 Chapter 5 Dynamic Cancer Model for Tumor Relapses et al. (1981), Lawless (1987), and Aalen and Husebye (1991). Second, the therapy does not have an effect on the patient, so the state of the patient is the same as just before the relapse. This is known in reliability literature as “minimal repair” and it has been studied Prentice et al. (1981), Brown and Proschan (1983), and Lawless (1987). We realize that some of the models mentioned above can be formulated using two different expressions for effective age function. As an example, some of the conditional models examined Prentice et al. (1981) (also called PWP in the introduction) can be formulated by organizing the data in calendar time (PWP-TT) (i.e., total time risk set or time from each unit’s entry into the observation set) or interoccurrence/gap time (PWP-GT) (i.e., gap time risk set or time since the previous event). Thus, the hazard function only differs in the at-risk process formulation. In the case of PWP-TT the effective age corresponds to Ei(s) = sand in the PWP-GP formulation to Ei(s) = s−SiN† i(s−). The choice between PWP-TT or PWP-GT depends on whether we are interested in the time that has elapsed since a patient entered the study or since the last relapse. Models which employ calendar time formulation assume that all interventions produce a minimal or no improvement in the patient. In medical terms, the disease is continuing in a stable manner. Models based on gap-time formulation assume that all interventions lead to perfect recovery for the patient, e.g. disease disappears, which is known as complete remission in the cancer literature. Other models that have been used in biomedical problems where neither ρ[·;α] nor Zare both identically unity are the following. Gail et al. (1980) cancer occurrence model is obtained if we take, Z= 1, ρ[N† i(s−); α] = max{α−N† i(s−),0}, where αis some real number, and λ0(s) = λ0, where λ0is some positive constant. In this model αcan be interpreted as an initial measure of the patient’s susceptibility to events, which is becoming weaker as the relapses accumulate. Shared frailty model (see Oakes, 1991) arises from Pe˜na and Hollander’s model by taking ρ[N† i(s−); α] = 1, Ei(s) = s−SN† i(s−), and putting some parametric distribution to the frailty component Z, such as a gamma or a lognormal distribution. In all examples previously mentioned, the expressions of Ei(s) do not depend on the number of relapses. In other words, the effect of treatment after each occurrence is always the same. However, we know that the effect of intervention after each relapse is not always the same. So, we need to define how to incorporate the effect of intervention upon relapses via a more general effective age. To do so, response to therapy will become crucial. The next section deals with this aspect. 5.2 Effective Age Process for Cancer Data 101 5.2 Effective Age Process for Cancer Data In Section 1.3.1 some examples about effective age in biomedical settings were discussed. Now, we will focus on cancer problem. Patients with indolent lymphomas, and in general with cancer, are monitored from the date of diagnosis to the time of death or loss to follow-up. At any evaluation, complete remission (CR) is defined as the disappearance of tumor masses and disease-related symptoms, as well as the normalization of the previous test and/or biopsies, lasting for at least one month. Partial remission (PR) is said to occur when measurable lesions have decreased by at least 50%. Patients not included in these categories are called non-responders (NR) (see Cheson et al., 1999). It is well-known that the response to the treatment is related to the time until the next relapse (see for instance Montoto et al., 2002, Weisdorf et al., 1992, or Davidge-Pitts et al., 1996), so it is reasonable to have a model which incorporates this information. However, neither AG, PWP, WLW, GSB, nor frailty models, which are mostly used in biomedical settings, have incorporated the effect of performed interventions upon event reoccurrences, though some of these methods could accommodate such information through the use of time-dependent covariates as will be illustrated in the example presented in Section 5.4. We propose using the response to therapy, defined as CR, PR, or NR, to define a model for the effective age for relapsing patients as follows. Consider a single patient and let {Aj:j= 0,1,2, . . .}be a sequence satisfying A0= 0, Aj=Aj−1+Ãj Y k=1 [1 −ψk]!Tj, j ≥1,(5.1) where ψ= (ψ1, ψ2, . . . , ψK)0,with ψj∈ {0, .5,1}and with the interpretation that ψj= 0 means that an NR (non-response) has occurred after the jth relapse, ψj= 1 means that a CR (perfect intervention) has occurred, while ψj=.5 means that a PR (partial remission) has transpired. The values of the ψjs can be assessed by the clinician(s) monitoring the patient. Our proposed effective age process for cancer relapse is E(s) = AN†(s−)+³s−SN†(s−)´.(5.2) The effective age (5.2) is a particular case of Kijima’s 1989 model II (Section 1.3.2), where in his model the ψjs are assumed to take any values in [0,1], whereas in our model we assume that they only take three possible values. Kijima’s aim was to model the situation where after 102 Chapter 5 Dynamic Cancer Model for Tumor Relapses each failure some repair is performed and the effectiveness of this repair could be quantified by a number that is between zero and one. We note that in cancer problems, we may also assess this “degree” of response according to a number in [0,1]; however, in realistic settings clinicians only need to know if a CR, PR, or NR was achieved to make a good therapy determination, hence our restriction of the possible values of the ψjs to the set {0, .5,1}. Dorado et al. (1997) also studied effective age functions that encompass the form (5.2). Note that if all responses are CR, i.e. ψi= 1, i = 1,2, . . ., then the effective age corresponds to gap time formulation, E(s) = s−SN†(s−)since all Ajs in (5.1) become 0. Similarly, if all responses are NR, the effective age corresponds to a calendar time formulation, E(s) = s. To demonstrate the notion of an effective age in biomedical settings, Figure 5.1 shows the effective age for a patient in a cancer study. This process between 0 (or S0) and S2corresponds to E(s) = s(calendar or elapsed time formulation). At the first event S1, treatment or intervention did not improve the disease status. In medical parlance, the patient did not respond to treatment, i.e., NR is achieved. After the second event, which occurred at S2, the patient responds perfectly to treatment, achieving a complete resolution of all clinical manifestation of the disease. It is considered a CR. In this case the effective age corresponds to E(s) = s−SN†(s−)(backward recurrence time). However, after the third event at time S3, the patient reverts to a state between a CR and a NR, that is, the patient experiences a little improvement or a PR to treatment. Finally, a progressive disease is observed for the fourth relapse at S4, possibly due to some complications. We have decided to show a progressive possibility despite we do not include it in our analysis because physicians may observe this types of responses. Thus, we indicate how this information should be incorporated in the model. Finally, the fifth failure which would have happened at S5 is not observed since the end of observational period τfor this hypothetical patient is less than S5. Consequently, the gap time for the fifth event is right-censored by τ−S4. . Next, we illustrate the effective age with a numerical example. Let us suppose that we observe a patient who receives an initial treatment at time 0, and this patient relapses four times at calendar times 30, 55, 100, and 150, and gets censored at calendar time 175. Thus, the gap times are 30, 25, 45, and 50, respectively. If we assume perfect repair model, e.g. the patient achieves CR after each intervention, the effective age at each recurrence will be 0, 30, 25, 45, and 50. On the other hand, if we assume minimal repair model, e.g. treatment is not effective at each recurrence or NR after each treatment, the effective age at each relapse will be 0, 30, 55, 100, and 5.2 Effective Age Process for Cancer Data 103 Figure 5.1: Pictorial representation of effective age vs. calendar time for a hypothetical unit in cancer settings. 150 (the same as calendar times). But the effect of the intervention at each recurrence need not be always the same. For example, let us suppose that the treatment after the first relapse does not improve the patient’s health at all, so we observe an NR. Then after the second relapse the patient achieves a CR. A PR is observed after third recurrence, and finally an NR is achieved at the fourth occurrence. Thus, if the patient has no improvement after the first intervention, the effective age will be the same as if the patient had a minimal repair. After the second relapse the patient has a perfect intervention. In that case the effective age corresponds to gap time or backward recurrence time. After the third recurrence the patient acquires some, but not total, improvement, so the effective age will be between those observed in the perfect and the minimal repair situations, say halfway. Finally, in the last relapse patient does not get better, so the effective age will start at a higher value and proceed possibly in a linear fashion. Under this hypothetical situation, the effective age would take the values 0, 30, 55, 45, and 72.5, with the effective age at the time of censoring (calendar time 175) being higher than 97.5. We will refer to Pe˜na and Hollander model with effective age in 5.1 as a dynamic cancer model. Because this model is just a special case of the general recurrent event model considered in Chapter 4, the procedures for estimating the parameters of this general class of models therefore applies to this cancer model. 104 Chapter 5 Dynamic Cancer Model for Tumor Relapses 5.3 Simulation Study 5.3.1 Simulation Design We have carried out simulation studies to examine empirically the properties of the parameter estimators described in Sections 5.1 and 5.2. In particular we study in the simulation: •the effect of sample size (n) •bias and variance of the estimators •the performance of the estimator of the baseline survival function ¯ F0(t) = exp{−Rt 0λ0(w)dw}, in terms of its bias function and root-mean-square-error (RMSE) function at specified time points •the consequences of miss-specifying the effective age function In the last item we consider the situation where the data have been generated by the model allowing for different responses after each intervention, but with the resulting data analyzed by assuming that the patients always achieve the same response (that is, either always CR or always NR) after each tumor reoccurrence. We mimicked the simulation study performed by Pe˜na et al. (2003). We point out that because the simulation was meant to cover general biomedical settings, some of the parameter values we considered may not be realistic for cancer settings, for instance, the mean number of events per patient, and effective age function. For the simulation, we considered an effective age function corresponding to the cancer model. That is, a patient can achieve a complete, a partial, or a null response depending on the vector ψ. We have assumed three different scenarios according to the following probability functions for ψ which takes values in the set {1=CR, 0.5=PR, 0=NR}:{(.8, .1, .1),(.3, .5, .2),(.1, .2, .7)}. Thus, in the first case, we assume that patients achieve CR with a probability of 80%, and PR or NR 10% of the time, respectively. These three sets of distributions allow us to cover three different scenarios: the first assumes that in a large majority of cases, perfect response is achieved after each relapse; the third has minimal response predominating; and the second distribution is an in-between scenario. For notation in the sequel, when we write p(ψ) = (p1, p2, p3) to indicate that 5.3 Simulation Study 105 the ψvalues are chosen such that 1 (= CR) occurs with probability p1, .5 (= PR) occurs with probability p2, and 0 (= NR) occurs with probability p3. To examine the impact of sample size, we select two values of n:n∈ {30,50}. These values are sufficient to study the limitations of the model in cancer settings since it is usual to have at least these sample sizes in this context. The censoring variables τi,i= 1, . . . , n, are generated according to same scheme as in previous simulation study. That is, a uniform distribution over [0, B], where Bis chosen such that under the assumption that patients always achieve a complete remission after relapses (i.e., E(s) = s−SN†(s−)) and with no accumulating events effect (i.e., α= 1). On average, there are approximately 5 events per patient since in cancer problems, it is difficult to find situations with more than this number of reoccurrences. For the baseline hazard function λ0, we choose a Weibull distribution, with unit scale parameter and shape parameter, γ, taking values in {.9,2}. Thus, we are able to study two different situations: one for an increasing baseline hazard rate function and the other one for a decreasing case. The impact of the accumulating number of relapses is assumed to be of form ρ(k, α) = αkform. We have selected α∈ {0.9,1,1.05}. Thus, we are able to study the case where an increasing number of relapses increases the time of the next relapse (that is, beneficial effect) which is the case when α=.9, as well as the case where there is no effect which is for α= 1, and the case where there is an adverse effect which is when α= 1.05, respectively. In order to take into account the effect of covariates, we have simulated a two-dimensional covariate vector (X1, X2). Then, to have both categorical and continuous covariates, X1has been simulated to have a Bernoulli distribution with success probability of .5 and X2was set to have a standard normal distribution. These covariates were generated to be stochastically independent. The regression coefficient vector (β1, β2) was set to (1,−1). Finally, the frailty component was generated under a gamma distribution with unit mean and variance 1/ξ. The parameter ξtook values in {2,6,∞}, with ∞corresponding to the absence of frailties. We performed 1,000 replications for each combination of simulation parameters. To create the bias and RMSE curves of the estimator of the baseline survivor function, we chose the time values that corresponded to the [0 : (.01) : .99] quantiles of the true baseline distribution function. In order to study the miss-specification of effective age function, we have also estimated the parameters for the dynamic cancer model assuming a model with an effective age function that considers that patients always achieve a perfect response (CR) and another where patients 106 Chapter 5 Dynamic Cancer Model for Tumor Relapses always achieve a minimal response (NR). Through these simulations we are able to highlight the importance of the effective age function in relation to either under- or over-estimation of the baseline survivor function. Given that we wanted to compare three different effective age functions (minimal, perfect, and cancer model) and given the three discrete distributions that specified the effective age process, two sample sizes, two hazard shapes, tree levels of event dependence, and three degrees of correlation, we conducted a total of 324 simulation experiments. 5.3.2 Simulation Results In the discussion of the simulation results, we will focus on the consequences of analyzing data from the cancer model when analyzed using models which always assume the same response, either always perfect or always minimal. Regarding distributional properties of the estimators of α,β, and η≡ξ/(1 + ξ), and the estimator of the baseline survivor function when the correct model is utilized, we only focus on the results presented in the tables and figures. Further details of this simulation are described in Pe˜na et al. (2003). Results of the simulations are shown in the following tables. Table 5.1 summarizes the mean values and standard deviations of the sampling distributions of the estimators of α,β1,β2, and ηfor α,n, and ξvarying in the sets mentioned above and for p(ψ) = (0.8,0.1,0.1). Figure 5.2 shows plots of the bias and RMSE curves for the non-parametric estimator of ¯ F0, where α= 0.9, p(ψ) = (0.8,0.1,0.1), and n∈ {30,50}. The lines of each plot represent the three different values of ξ: 2, 6, and ∞. Each plot frame contains the figure for two Weibull shape parameters, γ= 0.9, and γ= 2.0. Tables 5.2 and 5.3 present the summary of simulation results belonging to the effective age mis-specification analysis. Figure 5.3 shows the estimated baseline survivor function, the bias, and the RMSE curves calculated under effective age miss-specification, showing the effect of different effective age function chosen (always minimal or perfect response), the impact of different ξvalues, and for the three different cancer models analyzed. The results are for α= 0.9 and n= 30. As one may expect, when there is no miss-specification and when the sample size increases, the performance of the estimators of the finite-dimensional parameters and the baseline survivor function improved, as can be seen by noting that both bias and standard error decrease. We also notice that when the sample size is small, there is considerable over-estimation of η(Table 5.1). Examining the curves in Figure 5.2 we observe that the estimator of baseline survivor function is 5.3 Simulation Study 107 α γ ξ η n NC ˆµEv ˆµˆαˆσˆαˆµˆ β1ˆσˆ β1ˆµˆ β2ˆσˆ β2ˆη 0.9 0.9 2 0.67 30 0 2.87 0.896 0.035 1.033 0.415 −1.012 0.241 0.714 0.9 0.9 2 0.67 50 0 3.98 0.896 0.022 1.021 0.325 −1.024 0.174 0.69 0.9 0.9 6 0.86 30 0 6.13 0.895 0.032 1.02 0.327 −1.027 0.198 0.895 0.9 0.9 6 0.86 50 0 4.5 0.897 0.022 1.005 0.23 −1.015 0.141 0.877 0.9 0.9∞30 0 5.4 0.894 0.026 1.022 0.237 −1.022 0.152 0.9 0.9∞50 0 4.1 0.896 0.019 1.03 0.169 −1.022 0.1 0.9 2 2 0.67 30 0 7.13 0.905 0.019 0.981 0.302 −0.994 0.168 0.763 0.9 2 2 0.67 50 0 8.32 0.904 0.014 1.014 0.222 −0.985 0.121 0.731 0.9 2 6 0.86 30 0 6.8 0.903 0.017 0.983 0.206 −0.987 0.127 0.896 0.9 2 6 0.86 50 0 7.9 0.903 0.013 0.997 0.154 −1.001 0.09 0.883 0.9 2 ∞30 0 7.5 0.899 0.017 1.036 0.161 −1.02 0.099 0.9 2 ∞50 0 8.2 0.898 0.012 1.014 0.113 −1.016 0.067 1 0.9 2 0.67 30 10 2.87 0.994 0.032 1.003 0.444 −1.001 0.251 0.728 1 0.9 2 0.67 50 1 4 0.999 0.015 1.006 0.321 −0.99 0.192 0.713 1 0.9 6 0.86 30 5 5.2 0.995 0.029 1.033 0.353 −1.019 0.199 0.903 1 0.9 6 0.86 50 0 4.72 0.998 0.014 1.015 0.264 −1.009 0.155 0.884 1 0.9∞30 0 2.3 0.995 0.024 1.006 0.254 −1.031 0.165 1 0.9∞50 0 4.8 0.997 0.014 1.025 0.198 −1.019 0.121 1 2 2 0.67 30 0 5.07 1.009 0.027 0.988 0.327 −0.973 0.183 0.744 1 2 2 0.67 50 0 4.86 1.008 0.019 0.99 0.25 −0.981 0.137 0.73 1 2 6 0.86 30 1 4.97 1.006 0.026 1 0.262 −0.988 0.144 0.889 1 2 6 0.86 50 0 4.06 1.007 0.017 0.993 0.19 −0.991 0.11 0.867 1 2 ∞30 0 7.23 0.997 0.023 1.031 0.202 −1.023 0.118 1 2 ∞50 0 3.92 0.999 0.016 1.024 0.149 −1.014 0.094 1.05 0.9 2 0.67 30 5 6.53 1.05 0.017 1.021 0.421 −0.997 0.239 0.742 1.05 0.9 2 0.67 50 0 6.42 1.051 0.008 0.984 0.315 −0.991 0.184 0.707 1.05 0.9 6 0.86 30 6 2.03 1.051 0.016 0.989 0.317 −1.003 0.205 0.894 1.05 0.9 6 0.86 50 3 4.96 1.051 0.008 1.002 0.231 −0.994 0.154 0.889 1.05 0.9∞30 0 6.2 1.05 0.014 1.029 0.229 −1.036 0.155 1.05 0.9∞50 0 7.34 1.051 0.008 1.022 0.166 −1.01 0.106 1.05 2 2 0.67 30 1 5.77 1.055 0.022 0.988 0.334 −0.98 0.187 0.758 1.05 2 2 0.67 50 0 7.26 1.051 0.012 1.01 0.237 −1.006 0.138 0.72 1.05 2 6 0.86 30 0 8.63 1.053 0.02 1.008 0.246 −1.011 0.145 0.889 1.05 2 6 0.86 50 0 7.54 1.052 0.013 1.005 0.192 −1.008 0.105 0.872 1.05 2 ∞30 0 6.57 1.049 0.018 1.034 0.187 −1.012 0.119 1.05 2 ∞50 0 7.2 1.048 0.013 1.03 0.141 −1.023 0.082 Table 5.1: Summary of simulated means and standard deviations of the estimators of α,β, and η= ξ/(ξ+ 1). The true value of βis (1,−1). Results correspond to the case of p(ψ) = (0.8,0.1,0.1), and 1000 replications were done for each parameter combination. The others columns are: γ Weibull shape parameter, nsample size, NC number of replicates in which there was no model convergence; ˆµEv mean number per patient in all the simulation replications. 114 Chapter 5 Dynamic Cancer Model for Tumor Relapses MinimalaPerfectbCancerc α0.72 (.37) 0.90 (.15) 0.68 (0.21) Frailty ξ2.24 ∞1.36 ν0.45 8.97×10−80.73 Lesions Single 1 1 1 Localized 2.48 (0.93-6.60) 2.42 (0.90-6.50) 2.59 (0.88-7.62) >1 nodal site 4.47 (1.64-12.15) 3.26 (1.20-8.87) 4.55 (1.22-16.93) Generalized 3.24 (0.70-14.95) 2.69 (0.59-12.31) 3.09 (1.01-9.37) aEffective Age is E(s) = s. bEffective Age is E(s) = s−SN†(s−). cEffective Age is AN†(s−)+¡s−SN†(s−)¢. Table 5.4: Hazard ratios and confidence intervals at 95% (in parenthesis) for the probability of relapse depending on lesions involved at diagnosis for the non-Hodgkin’s lymphoma data set. Estimates obtained from the general model using three different effective ages processes. for further details). After preparing the data and including the treatment response as a dummy variable, the HR for variable lesions at diagnosis are: 2.46 (CI95% 1.06 to 5.77), 3.25 (CI95% 1.25 to 8.47), and 2.77 (CI95% 0.83 to 9.22), respectively. These results are similar to those obtained using the perfect repair model. Considering that the cancer model reveals that a frailty component is important, perhaps the AG model is not adequate since this model assumes that there is no heterogeneity among patients. Finally, we compare our results to those obtained using only time to first relapse in a Cox model. In that case, the HR are: 1.40 (CI95% 0.48 to 4.03), 2.76 (CI95% 0.95 to 8.03), and 3.24 (CI95% 0.28 to 37.50). Here, we ignore the information in the subsequent relapse times and we observe that this fact substantially affects the estimate of the coefficients and their statistical significance, especially in patients with localized lesions. The heterogeneity of the risk of relapse may be explained by subject-specific factors other than lesions involved at diagnosis, such as gender or delay between first treatment and first symptom. Thus, we now include all three covariates and compare the estimates of the regression coefficients from the cancer model with those obtained using some of the currently-used models (AG, AG with time-dependent covariates, WLW, and shared gamma frailty model). Table 5.5 shows these resulting HRs for the PCMZCL data. After adjusting for gender and delay, the variance of frailty decreased to 0.11 (1/8.85) indicating that the frailty is not necessary (likelihood ratio test 2(180.07 −179.40) = 1.34, p= 0.2476). Similarly, although the estimate of αdiffers from 1, it is not statistically significant (based on assumed asymptotic normality), so it seems that the 5.4 The non-Hodgkin’s lymphoma study 115 prior number of event occurrences does not have an impact. The results also indicate that shared frailty model gives lower risk estimates than those obtained using cancer model, while WLW method gives higher risk estimates than cancer model. Only models based on AG approach show similar results to those obtained using cancer model. In that case gender differences is statistically significant only if the AG model is chosen, while cancer model is the only one that indicates an statistically significant increased risk for those patients with generalized lesions as compared to patients with single lesions. The estimates of the cumulative hazard functions for the multiple event data for cancer model and for AG model with treatment response as a time-dependent covariate are shown in Figure 5.5. The solid lines correspond to hazard of relapse for patients with single lesions at diagnosis, obtained via ˆ Λ0(s) exp(ˆ β2¯ X2) and ˆ Λ0(s) exp(ˆ β1+ˆ β2¯ X2) for males and females, respectively, where ˆ Λ0(·) is the estimated cumulative baseline hazard. The dotted lines in this figure are for patients with generalized lesions which correspond to ˆ Λ0(s) exp(ˆ β2¯ X2+ˆ β5) and ˆ Λ0(s) exp(ˆ β1+ˆ β2¯ X2+ˆ β5) for males and females, respectively. The observed means are ¯ X2= 2.4 for males, and ¯ X2= 2.7 for females. These plots indicate that different risks are associated with the number of lesions involved at diagnosis, as is clear from their associated HRs. The AG and cancer model estimates of the hazard rate functions differ more for patients with generalized lesions than with single lesions. 116 Chapter 5 Dynamic Cancer Model for Tumor Relapses AG AG2 Shared Covariate Frailty Frailty WLW Frailty Cancera α- - - - .88 (.40) Frailty ξ24.51 ∞-∞8.85 ν0.04 5 ×10−7- 5 ×10−70.11 Gender Males 1 1 1 1 1 Females 2.01 2.01 1.83 1.73 1.84 (1.01-4.00) (1.02-3.92) (0.81-4.14) (0.88-3.40) (0.82-4.10) delay in years 0.99 1.01 1.02 1.04 0.99 (0.89-1.12) (0.89-1.12) (0.88-1.18) (0.77-1.40) (0.80-1.23) Lesions Single 1 1 1 1 1 Localized 3.70 3.83 5.23 3.24 3.57 (1.18-1.16) (1.23-11.9) (1.71-15.96) (1.05-9.96) (1.17-10.89) >1 nodal site 4.71 4.77 6.45 3.99 4.67 (1.62-13.7) (1.62-14.0) (2.37-17.56) (1.41-11.28) (1.25-17.4) Generalized 4.75 4.60 23.16 3.44 4.60 (0.92-24.4) (0.86-21.5) (5.02-106.9) (0.69-16.97) (1.30-16.3) aEffective Age is AN†(s−)+¡s−SN†(s−)¢. Table 5.5: Hazard ratios and confidence intervals at 95% (in parenthesis) for the probability of relapse for the PCMZCL data set. Estimates from the Andersen-Gill (AG), and Andersen-Gill with response to treatment after relapse as time-dependent covariate (AG2) including a frailty term, together with the estimates obtained form Wei, Lin and Weissfeld (WLW), Shared Gamma Frailty model, and dynamic cancer model. 5.4 The non-Hodgkin’s lymphoma study 117 0 20 40 60 80 100 120 0.0 0.5 1.0 1.5 2.0 Time (months) hazard Effective age function Cancer model AG model with response as time−dependent covariate Males, single lesions vs. generalized 0 20 40 60 80 100 120 0.0 0.5 1.0 1.5 2.0 Time (months) hazard Effective age function Cancer model AG model with response as time−dependent covariate Females, single lesions vs. generalized Figure 5.5: Estimates of cumulative hazard function for multiple events for non-Hodgkin’s lymphoma data set by sex and lesions involved at diagnosis, all cases evaluated at the mean value of delay between first treatment and first symptom. The blue line shows the hazard assuming the cancer model and red lines correspond to AG model which includes the response to treatment as a time-dependent covariate. Solid lines are patients with single lesions and dotted lines are for patients with generalized lesions at diagnosis. 118 Chapter 5 Dynamic Cancer Model for Tumor Relapses 5.5 R instructions for gcmrec package Parameter estimates for the dynamic cancer model can be obtained using gcmrec function and indicating where it is the information about effective age. In our case cancer=lymphoma$effage. > mod.can<-gcmrec(Survr(id,time,event)~as.factor(distrib), + data=lymphoma, s=1000, Frailty=TRUE , se="Jacknife", + cancer=lymphoma$effage) > mod.can Call: gcmrec(formula = Survr(id, time, event) ~ as.factor(distrib), data = lymphoma, s = 1000, Frailty = TRUE, se = "Jacknife", cancer = lymphoma$effage) coef exp(coef) se(coef) Jacknife z p as.factor(distrib)1 0.953 2.59 0.556 1.7146 0.086 as.factor(distrib)2 1.516 4.55 0.667 2.2738 0.023 as.factor(distrib)3 1.129 3.09 0.569 1.9842 0.047 General class model parameter estimates rho function: Alpha to k alpha (s.e. Jacknife): 0.683 (0.214) Frailty parameter, Xi (s.e. Jacknife): 1.36 Marginal log-likelihood= -176.27 n= 63 n times= 112 number of iterations: 74 EM steps Chapter 6 Concluding Remarks and Future Research 6.1 Conclusions The general class of model proposed by Pe˜na and Hollander (2004) have been demonstrated to be very useful to deal with recurrent event data. In particular, the dynamic cancer model developed in this PhD thesis has demonstrate to be very useful in analyzing indolent diseases in which a relapsing pattern is observed. We have to mention that although some physicians are interested in analyzing this type of data (see MacLaughlin, 2002), these models have received no attention, so far. The model we propose for analyzing cancer data, which includes as special cases many well-known models in survival analysis, is also important because it takes into account the effect of interventions which are performed after each event occurrence through the notion of an effective age, the possible weakening (or strengthening) effect of accumulating event occurrences, the possible presence of unobserved frailties that could be inducing correlations among the interevent times per unit, and the effect of observable covariates. Regarding procedures for estimating the parameters of this model, three different approaches have been described. One of them is based on the EM algorithm, while the other two used penalized likelihood inference. Regarding the general behavior of the model proposed, our simulations suggested that an under-specification of the model, in the sense of analyzing a data generated from the model with frailties using procedures developed from the model without frailties, could 119 120 Chapter 6 Concluding Remarks and Future Research have unacceptable consequences in that the resulting estimators will have non-negligible systematic biases. On the other hand, it was found that over-specification of the model may provide a robust method of analysis with an acceptable loss in efficiency. The application of the procedures to the bladder cancer data set highlights the importance of monitoring the effective age process. The simulations for the miss-specification of effective age by not incorporating information about the intervention effect in cancer settings, but instead assuming that always a complete or null response is achieved, has undesirable consequences because the resulting estimators of the finitedimensional parameters and the baseline survivor function are highly biased. Our application of the cancer model to an indolent lymphoma data set also highlights the need to incorporate information about the effect of intervention after each relapse. The main advantage of using the general class of models proposed by Pe˜na and Hollander (2004) with respect to analyze recurrent event data using other existing models, is the ability of incorporating information about the performed interventions after reoccurrences. In particular, we have shown that this model may be useful to analyze data set arising form cancer settings, in which the response to the treatment after each relapse is an important factor to predict new relapses. The dynamic cancer model may be used in a variety of applications when information about the response to intervention upon relapses can be obtained. The bladder cancer data set is an example. This data set, however, does not contain information about the effective age function and this leads to be the main limitation of this model. One possibility to solve this drawback is to use simple forms of the effective age, such as perfect or minimal repair formulations as it has been illustrated through this thesis. In this sense, the cancer model tries to model more complicated forms for the effective age function in which the response to the treatment after relapses are included in the model acting in the baseline hazard function. On the other hand, this PhD thesis has also addressed the problem of how to calculate confidence intervals for median survival time. We have proposed two different methods. One of them is based on asymptotic variances in the case that interoccurrence times are i.i.d. and another used bootstrap techniques. We have also studied several bootstrapping schemes to estimate the sampling distribution of median survival time estimators in the presence of recurrent event data and in consideration of the sum-quota data accrual which induces informative stopping and censoring. We proposed several resampling plans under the i.i.d. model and a correlated interoccurrence times model. From the simulations studies we carried out, we may conclude that 6.2 Future Research 121 the best bootstrapping scheme to estimate the median survival sample distribution under an i.i.d. model are just bootstrapping from observed data and non-parametric bootstraps (plans I, and II or III, respectively). For a correlated interoccurrence times (under a gamma frailty model), both semiparametric plans (VI and VII) are the best ones. Plan IV, which is anchored in using the WC (1999) estimator of the inter-event survivor function appears to offer a robust procedure when uncertain about the model that generated the data. Based on the simulation studies, it appears that bootstrapping from the empirical distribution of the monitoring times do not provide improvements. 6.2 Future Research Maybe the main important future work regarding confidence interval for median survival time is to prove asymptotic convergence of bootstrap procedures. There are still many important questions that need to be examined with regards to the general model proposed by Pe˜na and Hollander (2004). The first is the ascertainment of asymptotic properties of the estimators, such as their asymptotic normality or the weak convergence to a Gaussian process of a properly normed estimator of the baseline survivor function. The resolution of this asymptotic problem may require methods utilized in Murphy (1994, 1995) and Parner (1998). Some asymptotic results for the general class of models when baseline hazard function is parametrically specified can be found in Kvam and Pe˜na (2003). Through such asymptotic analysis we will be able to obtain expressions for approximating analytically the standard errors of the estimators and compare them with those obtained using penalized approaches. Another issue of importance is whether bootstrapping methods could be utilized to obtain standard errors of estimators for the purpose of constructing confidence intervals and/or bands for the parameters. The problem of how to validate this class of models after it has been fitted to a specific data set is another open problem, and calls for suitable goodness-of-fit and model validation procedures. Regarding the dynamic cancer model, new applications for the cancer model may require modifications of our formulation. One important consideration is the value of A0in the effective age (see equation (5.2)). Because all our lymphoma patients achieved CR at first treatment, we were able to use this date as study origin, assured that all patients had the same initial status with A0= 0. In other situations, however, CR may not be achievable at first treatment for all 122 Chapter 6 Concluding Remarks and Future Research subjects. In this case, A0is not uniformly zero for all subjects, and disease status after first treatment may be assessed to assign a positive value for A0. Another aspect of our work that requires consideration in future research is that the time between treatment and assessment of the response to treatment (the ψjs) following each relapse may not be negligible, contrary to our context and earlier developments in reliability. The delay between application of treatment and evaluation of patient response is widely recognized. However, in cancer studies, at least in hematological diseases, this may not be a problem since all patients are routinely monitored for one month following administration of therapy. Bibliography O. Aalen. Heterogeneity in survival analysis. Stat. Med., 7:1121–1137, 1988. O. Aalen and E. Husebye. Statistical analysis of repeated events forming renewal processes. Stat. Med., 10:1227–1240, 1991. P. Andersen, O. Borgan, R. Gill, and N. Keiding. Statistical Models Based on Counting Processes. Springer-Verlag, New York, 1993. P. Andersen and R. Gill. Cox’s regression model for counting processes: a large sample study. Ann. Statist., 10:1100–1120, 1982. U. Barai and N. Teoh. Multiple statistics for multiple events, with application to repeated infectious in the growth factor studies. Stat. Med., 16:941–949, 1997. M. Barcel´o. Modelos marginales y condicionales en el an´alisis de supervivencia multivariante. Gac. Sanit., 16 (Supl. 2):59–68, 2002. R. Beran. Estimated sampling distributions: the bootstrap and competitors. Annals of Statistics, 10:212–225, 1982. P.J. Bickel and D.A. Frieedman. Some asymptotic theory for the bootstrap. Annals of Statistics, 9:1196–1217, 1981. H. Block, W. Borges, and T. Savits. Age-dependent minimal repair. J. Appl. Prob., 22:51–57, 1985. R. Brent. Algorithms for Minimization without Derivatives. Englewood Cliffs, NJ: Prentice-Hall, 1973. N. Breslow. ”discussion of professor cox’s paper”. J. Royal Statist. Soc., B, 34:216–217, 1972. 123