Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3045702

Formats

Article sections

- Abstract
- 1 Introduction
- 2 A Clinical Study
- 3 Notation and Model
- 4 A Phase III Study of Prostate Cancer
- 5 Results
- 6 Discussion
- Supplementary Material
- References

Authors

Related links

Biometrics. Author manuscript; available in PMC 2011 February 27.

Published in final edited form as:

Published online 2009 June 8. doi: 10.1111/j.1541-0420.2009.01276.x

PMCID: PMC3045702

NIHMSID: NIHMS117723

Song Zhang: ude.nretsewhtuostu@gnahZ.gnoS

The publisher's final edited version of this article is available at Biometrics

See other articles in PMC that cite the published article.

We consider inference for data from a clinical trial of treatments for metastatic prostate cancer. Patients joined the trial with diverse prior treatment histories. The resulting heterogenuous patient population gives rise to challenging statistical inference problems when trying to predict time to progression on different treatment arms. Inference is further complicated by the need to include a longitudinal marker as a covariate. To address these challenges, we develop a semi-parametric model for joint inference of longitudinal data and an event time. The proposed approach includes the possibility of cure for some patients. The event time distribution is based on a non-parametric Pólya tree prior. For the longitudinal data we assume a mixed effects model. Incorporating a regression on covariates in a non-parametric event time model in general, and for a Póolya tree model in particular, is a challenging problem. We exploit the fact that the covariate itself is a random variable. We achieve an implementation of the desired regression by factoring the joint model for the event time and the longitudinal outcome into a marginal model for the event time and a regression of the longitudinal outcomes on the event time, i.e., we implicitly model the desired regression by modeling the reverse conditional distribution.

We discuss inference for data from a phase III clinical trial on treatments of metastatic prostate cancer. The challenges include patient heterogeneity due to prior treatment history and the need to include a regression on prostate specific antigen (PSA) as a longitudinal marker. We constuct a semi-parametric Bayesian model to address these challenges. It implements joint inference on event time and longitudinal observations, with the possibility that some patients are cured.

Let *T* be the event time and *Y* be the longitudinal covariate. Most existing approaches are based on factoring the joint model as *P*(*T, Y*) = *P*(*Y*)*P*(*T* | *Y*). The first factor is the longitudinal submodel *P*(*Y*), typically assumed to be a mixed model. The second factor is the survival submodel *P*(*T* | *Y*). In the following discussion, we use the terms event time, survival time, time to progression and failure time exchangably. There is an extensive literature on the joint modeling of longtitudinal and event time data without a cured fraction (De Gruttola and Tu, 1994; Tsiatis et al., 1995; Lavalley and De Gruttola, 1996; Wulfsohn and Tsiatis, 1997; Dafni and Tsiatis, 1998; Henderson et al., 2000; Xu and Zeger, 2001; Lin et al., 2002; Ibrahim et al., 2004). A review can be found in Tsiatis and Davidian (2004). Less work has been published on the joint modeling of longitudinal and event time data with cure. Law et al. (2002) proposed a model with the longitudinal process described by an exponential-decay-exponential-growth model and a mixture model to accomodate cure. The imputed values of the longitudinal measurements are covariates in a proportional hazard model. Brown and Ibrahim (2003) and Chen et al. (2004) assume event times to arise from the development of an unobserved number of metastatis-competent tumer (MCT) cells, modeled by a Poisson distribution. Subjects with zero MCT cells constitute the cure group. Yu et al. (2004) provide a recent review of joint longitudinal-survival-cure models.

Specific to modeling PSA, Pauler and Finkelstein (2002) propose a joint analysis using a change-point regression model for PSA trajectory, and a Cox model for cancer recurrence time with time-dependent covariates including functions of longitudinal parameters and imputed PSA mean function. Lin et al. (2002) considered a latent class model to uncover subpopulation structure for both PSA trajectories and a survival outcome. Given latent class membership, the longitudinal marker and outcome are assumed independent. The model assumes class-specific baseline hazard functions and accommodates possibly time-dependent covariates. Yu et al. (2008) investigated individual prediction in prostate cancer studies using a joint longitudinal-survival-cure model. A logistic model is specified for the probability of an individual being in the susceptible group and separate nonlinar mixed effect models are assumed for the cure and susceptible groups. The event-time process is modeled by a proportional hazard model with time-dependent covariates including the current slope and current value of the PSA trajectory.

In the joint analysis of longitudinal and event time data, most researchers assume parametric or semi-parametric models for *P*(*T* | *Y*). However, it is difficult to implement non-parametric models for *P*(*T* | *Y*) because most non-parametric models do not allow straightforward incorporation of a regression on covariates. We propose to use the alternative factorization, *P*(*T, Y*) = *P*(*T*)*P*(*Y* | *T*). We proceed under the Bayesian paradigm. Choosing a non-parametric model for *P*(*T*) is the traditional problem of non-parametric inference for an event time. The model *P*(*Y* | *T*) is part of a convenient factorization of the joint model. It should not be interpreted as a predicting model of *Y* by future outcome *T*. We propose a mixed effect model for *P*(*Y* | *T*), with *T* conveniently included as a univariate covariate. The opposite, including a high-dimensional covariate *Y* in a non-parametric model for *T*, would pose challenging technical problems (This is distinct from the related problem of non-parametrically modeling a mean function *E*(*T* | *Y*), e.g., kernel smoothers). Both factorizations lead to a joint model, *P*(*T, Y*), describing the dependence between *T* and *Y*. It is this joint model that ultimately allows improved prediction of the event time given repeated measurements of the marker. Under the factorization *P*(*T*) *P*(*Y* | *T*) the desired regression *P*(*T* | *Y*) is not explicitly parametrized, but implied by Bayes theorem. In a parametric model and using maximum likelihood estimation the same factorization is used in Pawitan and Self (1993) to model longitudinal markers for AIDS patients. They assume Weibull models for the infection time and disease occurrence and a generalized linear model for the longitudinal measurements of T4 counts and T4/T8 ratio, with the intercept and slope being functions of the event times.

We use a Polya tree (PT) prior to model the event time distribution. The reasons for this modeling choice are the possibility to model multimodal distributions reflecting the diversity of the patient population, the computational simplicity, and the easy a priori centering at a parametric model. A PT prior can be constructed to give probability one to the set of continuous or absolutely continuous probability measures (Lavine, 1992).

Muliere and Walker (1997) implemented PT models in a survival analysis. Walker and Mallick (1997; 1999) applied PT priors in hierarchical generalized linear models, frailty models, and accelerated failure time models. Hanson and Johnson (2002) developed a general approach to modeling residual distributions with a mixture of PT. Neath (2003) used PT to model censored data. Paddock et al. (2003) developed randomized PT models, which uses random partitions to smooth out the effect of partitions on posterior inference. Hanson et al. (2007) used mixtures of PTs to construct a joint model for time-dependent covariates and survival time. They introduced flexible PT priors for the baseline distributions in the Cox model, the proportional odds model, and an accelerated failure time model accommodating time-dependent covariates. Their approach uses the factorization *P*(*T, Y*) = *P*(*Y*)*P*(*T* | *Y*). See Hanson (2006) for a review of recent development in finite PT models.

Androgen ablation (AA) is the preferred treatment for metastatic prostate cancer. AA therapy alters the natural history of the disease by disrupting the growth promoting effects mediated by androgen receptor signaling, which is usually accomplished by medical suppression of testicular endocrine function. Unfortunately, most patients with clinically detectable metastatic disease when the AA therapy started will eventually progress to androgen independent prostate cancer (AIPC). AIPC is a relentlessly progressive disease state, and is the cause of death for the vast majority of men in whom it develops. By this mechanism, prostate cancer leads to an annual death toll of more than 27,000 men in the United States.

To date, no treatment has been found to be curative for AIPC, and it is only fairly recently that some therapies are shown to alter the natural history of the disease. A chemotherapy demonstrated a survival advantage over historical results in a phase II trial conducted at M.D. Anderson Cancer Center (Ellerhorst et al., 1997). This therapy, dubbed KA/VE, treats patients with ketoconazole and doxorubicin alternating with vinblastine and estramustine.

In this paper we analyze data from a phase III trial at M.D. Anderson Cancer Center that comapred conventional AA therapy versus AA therapy plus three 8-week cycles of KA/VE. The aim of this trial was to investigate whether better clinical benefit can be achieved by applying the chemotherapy “early”, i.e., before the metastatic prostate cancer develops into the far-advanced AIPC. The two treatment arms are denoted by AA and CH, respectively. The patient population includes metastatic prostate cancer patients whose high risk of developing AIPC justifies long-term, sustained, androgen ablation. The primary endpoint is the time to progression (TTP) to AIPC, which is diagnosed by the following criteria: 1) Symptoms attributed by the treating physician to reflect progressive cancer; 2) Radiographic progression; 3) Rising PSA, with value greater than 1 and doubling time < 9 months; 4) Treatment with chemotherapy. The first 3 also require demonstration of testosterone < 50 and withdrawal of antiandrogens. More details about the clinical trial can be found in Millikan et al. (2008).

Besides TTP, we also observed the longitudinal measurements of PSA level from each patient. Carter et al. (2006) demonstrated that PSA velocity is associated with prostate cancer death even 10–15 years before diagnosis. To further improve the understanding of this important marker we propose to build a joint model of TTP and PSA.

To statisticians, a challenge posed by this clinical trial is the considerable heterogeneity among patients. Before coming to M.D. Anderson Cancer Center these patients had been treated by different physicians with different therapies at different institutions. These differences might have a long-term impact on the development of prostate cancer. Second, there is no completely satisfactory way to define “early” in the natural history of metastatic prostate cancer. As a practical solution, the clock start of the trial is defined as the initiation of the AA therapy. Thus at the beginning of the trial, the true stage of cancer might not be exactly the same for each patient.

We use *v* = 1, 2 to denote the two treatment arms (1 for CH and 2 for AA). Let *n _{v}* be the number of subjects in each arm. For the

We define the sampling model for the observed data (* y_{vi}, t_{vi}, d_{vi}*) from each patient. If

If ω_{vi} = 0, the subject is at risk of developing AIPC. We assume *T _{vi}* to be a random sample from distribution

Given *T _{vi}*, the longitudinal measurements

$$\begin{array}{ccc}{L}_{\mathit{\text{vi}}1}\hfill & =[{\mathit{y}}_{\mathit{\text{vi}}}\phantom{\rule{thinmathspace}{0ex}}|\phantom{\rule{thinmathspace}{0ex}}{T}_{\mathit{\text{vi}}},\mathbf{\Psi}][{T}_{\mathit{\text{vi}}}={t}_{\mathit{\text{vi}}}\phantom{\rule{thinmathspace}{0ex}}|\phantom{\rule{thinmathspace}{0ex}}{\omega}_{\mathit{\text{vi}}}=0,{G}_{v}][{\omega}_{\mathit{\text{vi}}}=0\phantom{\rule{thinmathspace}{0ex}}|\phantom{\rule{thinmathspace}{0ex}}{p}_{v}],\hfill & \text{for}\phantom{\rule{thinmathspace}{0ex}}{d}_{\mathit{\text{vi}}}=1,\hfill \\ {L}_{\mathit{\text{vi}}0}\hfill & =[{\mathit{y}}_{\mathit{\text{vi}}}\phantom{\rule{thinmathspace}{0ex}}|\phantom{\rule{thinmathspace}{0ex}}{T}_{\mathit{\text{vi}}},\mathbf{\Psi}]I({T}_{\mathit{\text{vi}}}>{t}_{\mathit{\text{vi}}})[{T}_{\mathit{\text{vi}}}\phantom{\rule{thinmathspace}{0ex}}|\phantom{\rule{thinmathspace}{0ex}}{\omega}_{\mathit{\text{vi}}},{G}_{v}][{\omega}_{\mathit{\text{vi}}}\phantom{\rule{thinmathspace}{0ex}}|\phantom{\rule{thinmathspace}{0ex}}{p}_{v}],\hfill & \text{for}\phantom{\rule{thinmathspace}{0ex}}{d}_{\mathit{\text{vi}}}=0.\hfill \end{array}$$

(1)

Note that *L*_{vi0} is an augmented likelihood with latent variable *T _{vi}* and ω

For *L*_{vi0}, the two values taken by ω_{vi} lead to two models of different dimensions. If ω_{vi} = 0, *T _{vi}* is a random parameter with prior

We assume prior independence, $[\mathbf{\Psi},{p}_{v},{G}_{v},v=1,2]=[\mathbf{\Psi}]\xb7{\displaystyle {\prod}_{v=1}^{2}[{p}_{v}]\xb7[{G}_{v}]}$. The prior specification [**Ψ**] and posterior inference for mixed models of repeated measurements have been discussed extensively. See, for example, Ibrahim et al. (2004). For priors [*p _{v}*],

For the unknown survival distribution *G _{v}*, we consider two choices. The first choice is a parametric model assuming

A PT prior is indexed by two hyper-parameters: a nested sequence of partitions Π = {*B*_{0}, *B*_{1}, *B*_{00}, *B*_{01}, …, *B*_{ε0},*B*_{ε1}, …} of the sample space *S*, with *S* = *B*_{0} *B*_{1} and *B*_{ε} = *B*_{ε0} *B*_{ε1}; and parameters = {α_{0}, α_{1}, α_{00}, α_{01}, …, α_{ε0}, α_{ε1}, …}. Here ε = ε_{1} ε_{m} denotes a binary sequence of length *m*. We can center the PT prior around a given distribution , by setting α_{ε0} = α_{ε1} and defining *B*_{ε} at level *m* to coincide with quantiles ^{−1}(*k*/2^{m}), *k* = 0, 1, , 2^{m}. The parameter has a similar role as the precision parameter in a Dirichlet process prior. Berger and Guglielmi (2001) considered a family of the form α_{ε1, , εm} = *c*·ρ(*m*), where ρ(*m*) = *m*^{2}, *m*^{3}, 2^{m}, 4^{m} or 8^{m}, and *c* > 0 is a constant. In general, any ρ(*m*) such that $\sum}_{m=1}^{\mathrm{\infty}}{\rho (m)}^{-1}<\mathrm{\infty$ guarantees the PT to be absolutely continuous. See Hanson (2006) for a recent review. More details are presented in Web Appendix B.

To facilitate discussion, we define the following notation. The set of observed and unobserved TTPs under treatment *v* are denoted by ${\mathit{t}}_{v}^{1}=\{{t}_{\mathit{\text{vi}}}:{d}_{\mathit{\text{vi}}}=1\}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{\mathit{T}}_{v}^{0}=\{{T}_{\mathit{\text{vi}}}:{d}_{\mathit{\text{vi}}}=0\}$, respectively. We also define ${\mathit{\omega}}_{v}^{0}=\{{\omega}_{\mathit{\text{vi}}}:{d}_{\mathit{\text{vi}}}=0\}$ to be the set of unknown indicators of cure. Without loss of generality, we assume that *d _{vi}* = 0 for

$$[\mathbf{\Lambda}\phantom{\rule{thinmathspace}{0ex}}|\phantom{\rule{thinmathspace}{0ex}}\mathit{Y},\mathit{t},\mathit{d}]\mathrm{\propto}{\displaystyle \prod _{v=1}^{2}\left\{\right({\displaystyle \prod _{i=1}^{{n}_{v0}}{L}_{\mathit{\text{vi}}0}}{\displaystyle \prod _{i={n}_{v0}+1}^{{n}_{v}}{L}_{\mathit{\text{vi}}1}})\phantom{\rule{thinmathspace}{0ex}}[{p}_{v}][{G}_{v}]\}[\mathbf{\Psi}],}$$

(2)

where (* Y, t, d*) = {(

Before proceeding with posterior MCMC, we analytically marginalize (2) with respect to *G _{v}*. Recall that each subject with ω

We compare the proposed model with four natural alternatives. Details of the competing models and results are described later, in Section 5. We use the conditional predictive ordinates (CPO) proposed by Gelfand et al. (1992) to compare different models. The CPO for subject *i* in group *v* (henceforth subject (*v, i*)) is defined as the posterior predictive distribution evaluated for the observation from subject (*v, i*), conditional on all the data minus the response from subject (*v, i*). Formally, letting (**Y**_{(−vi)}, **t**_{(−vi)}, **d**_{(−vi)}) = (* Y, t, d*) \ (

We return to the clinical trial from Section 2. The phase III trial for advanced prostate cancer had a total enrollment of 286 patients, with *n*_{1} = 137 in the CH arm and *n*_{2} = 149 in the AA arm. Starting from the diagnosis of prostate cancer, the PSA level of each patient was monitored for up to 10 years. On average, about 30 PSA measurements were collected from each patient. We use *y _{vij}* (

The horizontal axis indicates years after treatment. The censoring times are marked by +. “K-M Est.” denotes Kaplan-Meier estimates. “Model Est.” denotes estimates based on model *M*_{1}, where TTP are assumed to arise from **...**

PSA level normally increases as the prostate enlarges with age. When prostate cancer develops, however, it increases much faster. The typical effect of a treatment on PSA level is a sharp drop in PSA level immediately after the treatment. Then gradually, the body adjusts to offset the treatment effect, and the PSA level bounces back. The speed of rebound depends on the progress of cancer. Web Figure 1 plots the longitudinal profiles of four randomly selected patients. Note the variability among the profiles. Exploratory analysis indicates a negative correlation between the PSA slope and TTP. Based on these considerations, the longitudinal submodel [* y_{vi}* |

$${f}_{\mathit{\text{vij}}}({s}_{\mathit{\text{vij}}})={\theta}_{0\mathit{\text{vi}}}+{\theta}_{1\mathit{\text{vi}}}{s}_{\mathit{\text{vij}}}+{\gamma}_{2\mathit{\text{vi}}}({e}^{-{\phi}_{0\mathit{\text{vi}}}{({s}_{\mathit{\text{vij}}}-{u}_{\mathit{\text{vi}}0})}^{+}}-1)+{\eta}_{v}({e}^{-{\phi}_{1v}{({s}_{\mathit{\text{vij}}}-{u}_{\mathit{\text{vi}}1})}^{+}}-1)+{\gamma}_{1v}{({s}_{\mathit{\text{vij}}}-{u}_{\mathit{\text{vi}}1})}^{+}+({\theta}_{1\mathit{\text{vi}}}+{\gamma}_{1v})({e}^{-{\xi}_{v}{T}_{\mathit{\text{vi}}}}-1){({s}_{\mathit{\text{vij}}}-{u}_{\mathit{\text{vi}}1})}^{+},$$

(3)

where (*x*)^{+} = *x* if *x* > 0, and (*x*)^{+} = 0 otherwise. We assume independent normal residuals, ${e}_{\mathit{\text{vij}}}\stackrel{\mathit{\text{iid}}}{\sim}N(0,{\sigma}^{2})$. The first two terms define a line with intercept θ_{0vi} and slope θ_{1vi}, describing the baseline linear trend of PSA over age. The coefficients are subject-specific. Parameters η_{v} and _{1v} model the size and the slope of the drop after the intervention with CH or AA. As age *s _{vij}* moves beyond

As for the PT priors, *G _{v}* ~

The mixture probability *p _{v}* is assumed to be

To validate the proposed model we consider comparisons with four alternative models. Let *M*_{1} denote the proposed model (2). The second model, *M*_{2}, is also based on the factorization *P*(*T, Y*) = *P*(*T*)*P*(*Y* | *T*), with *P*(*Y* | *T*) as in (3), but *P*(*T*) being fully parametric. We assume a Weibull regression model for (*T _{vi}* | ω

$${f}_{\mathit{\text{vi}}}({s}_{\mathit{\text{vij}}})={\theta}_{0\mathit{\text{vi}}}+{\theta}_{1\mathit{\text{vi}}}{s}_{\mathit{\text{vij}}}+{\gamma}_{2\mathit{\text{vi}}}\{{e}^{-{\phi}_{0\mathit{\text{vi}}}{({s}_{\mathit{\text{vij}}}-{u}_{\mathit{\text{vi}}0})}^{+}}-1\}+{\eta}_{v}\{{e}^{-{\phi}_{1v}{({s}_{\mathit{\text{vij}}}-{u}_{\mathit{\text{vi}}1})}^{+}-1}\}+{\gamma}_{1\mathit{\text{vi}}}{({s}_{\mathit{\text{vij}}}-{u}_{\mathit{\text{vi}}1})}^{+}$$

and *e _{vij}* ~

$${h}_{\mathit{\text{vi}}}(t)={h}_{v0}(t)\phantom{\rule{thinmathspace}{0ex}}\text{exp}[{\zeta}_{1v}{f}_{\mathit{\text{vi}}}({u}_{\mathit{\text{vi}}1}+t)+{\zeta}_{2v}{f}_{\mathit{\text{vi}}}^{\prime}({u}_{\mathit{\text{vi}}1}+t)],$$

(4)

where *h*_{v0}(*t*) is the baseline hazard and (ζ_{1v}, ζ_{2v}) are scaling parameters. Under *M*_{4}, we model *h*_{v0}(*t*) as a piecewise constant function of *J* = 8 steps. For 0 < *q*_{1} < *q*_{2} < < *q*_{J−1} < ∞, we assume *h*_{v0}(*t*) = κ_{v1} if *t* ≤ *q*_{1}, *h*_{v0}(*t*) = κ_{v2} if *q*_{1} < *t* ≤ *q*_{2}, , and *h*_{v0}(*t*) = κ_{vJ} if *t* > *q*_{J−1}. Gamma priors are assumed for κ_{vj}. More details can be found in Ibrahim et al. (2004). Under *M*_{5}, we model *h*_{v0}(*t*) by a Weibull hazard, with Gamma priors for the scale and shape parameters. The priors of the other parameters are specified as in *M*_{1}.

The estimated LPML under *M*_{1} through *M*_{5} are 4833.6, 5115.2, 5007.9, 4889.1, and 5155.0, respectively. Clearly *M*_{1} achieves the best performance. The nonparametric PT model allows the density function to deviate from the form imposed by the Weibull assumption. Assuming a cure group further improves the model fit. Model *M*_{4} has the second best performance, which indicates that the PSA trajectory does play an important role in prostate cancer progression. The inferior performance of *M*_{5} suggests that the Weibull hazard assumption might be too restrictive for our data. We further validate the survival and cure aspect of the model based on subject specific martingale residuals (Barlow and Prentice, 1988; Therneau et al., 1990; Lin et al., 2002). The residuals are scattered horizontally over age (with three outliers), suggesting no evidence against the proposed model. The residual plot is shown in Web Figure 2.

The estimated cure probabilities *p _{v}* (

Under model *M*_{1}, different PSA profiles lead to different posterior distributions of *T _{vi}*. In Figure 3 we compare for four patients with censored TTP the PSA profiles (1st column) and the estimated posterior probability of “cure”

Web Figure 1 plots the longitudinal PSA profiles of four patients together with fitted values. Table 1 lists the posterior means and standard deviations of some parameters in *M*_{1}. The posterior estimates of ξ_{v} (*v* = 1, 2) are practically identical, implying that the impact of TTP on the trajectory of PSA profiles are similar across the two treatments. The estimates of γ_{1v} indicate that the PSA profiles of patients in the AA arm on average have an increased slope after treatment. The level and slope of the drop in PSA after the CH/AA treatment are modeled by *l _{v}*(

Given a currently observed PSA profile, we can use the proposed method to obtain the predictive distribution of TTP, which provides a good assessment of progression risk. This predictive distribution can be continuously updated with additional PSA measurements. We demonstrate this learning process in Figure 4. The left panel plots the PSA profiles of two hypothetical patients from the AA arm. Each point denotes a PSA measurement. The two patients have their PSA level measured at the same time points. Within the first two years the two PSA profiles are identical, and then they deviate: the first patient’s PSA level stays low, while the second gradually rises. The center panel shows the continuously updated posterior estimates of *P*(ω = 1 | **y**_{−t}), with **y**_{−t} being the accumulated PSA measurements up to the time of assessment. We interpret *P*(ω = 1 | **y**_{−t}) as the individual probability of long term survival. The right panel shows the continuously updated posterior estimates of *E*(*T* | ω = 0, **y**_{−t}).

We conducted a sensitivity analysis to explore the impact of *t _{c}* and

The proposed model allows researchers to relax parametric assumptions on the survival submodel imposed by existing methods. An important limitation is that *P*(*T, Y*) = *P*(*T*)*P*(*Y* | *T*) does not explicitly state how *T* is affected by *Y*. Given a particular longitudinal profile, we need to carry out posterior simulation to learn about the posterior survival distribution given *Y*. The proposed approach can readily be generalized to problems with more than two treatments. The longitudinal data model (3) is appropriate for the discussed application to the prostate cancer trial. In general, any well specified model with a regression on the event time could be used.

We thank the two reviewers and associate editor for their constructive suggestion. We thank Randall Millikan for supplying the data set for analysis. Our work was partially supported by the NIH CTSA Grant UL1 RR024982, the Cancer Center Support Core Grant CA16672, SPORE in prostate cancer grant CA90270 from the National Cancer Institute, National Institute of Health.

**Supplementary Materials**

Web Appendices and Figures referenced in the paper are available under the Paper Information link at the Biometrics website http://www.biometrics.tibs.org.

- Barlow WE, Prentice RL. Residuals for relative risk regression. Biometrika. 1988;75:65–74.
- Berger JO, Guglielmi A. Bayesian and conditional frequentist testing of a parametric model versus nonparametric alternatives. Journal of the American Statistical Association. 2001;96(453):174–184.
- Brown ER, Ibrahim JG. Bayesian approaches to joint cure-rate and longitudinal models with applications to cancer vaccine trials. Biometrics. 2003;59(3):686–693. [PubMed]
- Carlin BP, Chib S. Bayesian model choice via Markov chain Monte Carlo methods. Journal of the Royal Statistical Society, Series B: Methodological. 1995;57:473–484.
- Carter B, Ferrucci L, Ketterman A, Landis P, Wright J, Epstein JI, Trock B, Metter J. Detection of life-threatening prostate cancer with prostate-specific antigen velocity during a window of curability. Journal of the National Cancer Institute. 2006;98:1521–1527. [PMC free article] [PubMed]
- Chen M-H, Ibrahim JG, Sinha D. A new joint model for longitudinal and survival data with a cure fraction. Journal of Multivariate Analysis. 2004;91(1):18–34.
- Dafni UG, Tsiatis AA. Evaluating surrogate markers of clinical outcome when measured with error. Biometrics. 1998;54:1445–1462. [PubMed]
- De Gruttola V, Tu X. Modelling progression of CD4-Lymphocyte count and its relationship to survival time. Biometrics. 1994;50:1003–1014. [PubMed]
- Dellaportas P, Forster JJ, Ntzoufras I. On bayesian model and variable selection using mcmc. Statistics and Computing. 2002;12(1):27–36.
- Ellerhorst J, Tu S, Amato R, Finn L, Millikan R, Pagliaro L, Jackson A, Logothetis C. Phase II trial of alternating weekly chemohormonal therapy for patients with androgen-independent prostate cancer. Clinical Cancer Research. 1997;3:2371–2376. [PubMed]
- Gelfand A, Dey D, Chang H. Model determination using predictive distribution with implementation via sampling-based methods (with discussion). Bayesian Statistics 4 – Proceedings of the Fourth Valencia International Meeting; Oxford University Press.1992.
- Green PJ. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika. 1995;82:711–732.
- Hanson T, Branscum A, Johnson W. Technical report. University of Minnesota; 2007. Joint modeling of longitudinal and survival data using mixtures of polya trees.
- Hanson TE. Inference for mixtures of finite Polya tree models. Journal of the American Statistical Association. 2006;101(476):1548–1565.
- Hanson T, Johnson WO. Modeling regression error with a mixture of Polya trees. Journal of the American Statistical Association. 2002;97(460):1020–1033.
- Henderson R, Diggle P, Dobson A. Joint modelling of longitudinal measurements and event time data. Biostatistics (Oxford) 2000;1(4):465–480. [PubMed]
- Ibrahim JG, Chen M-H, Sinha D. Bayesian methods for joint modeling of longitudinal and survival data with applications to cancer vaccine trials. Statistica Sinica. 2004;14(3):863–883.
- Lavalley MP, De Gruttola V. Models for empirical Bayes estimators of longitudinal CD4 counts (Disc: P2337–2340) Statistics in Medicine. 1996;15:2289–2305. [PubMed]
- Lavine M. Some aspects of Polya tree distributions for statistical modelling. The Annals of Statistics. 1992;20:1222–1235.
- Law NJ, Taylor JMG, Sandler H. The joint modeling of a longitudinal disease progression marker and the failure time process in the presence of cure. Biostatistics (Oxford) 2002;3(4):547–563. [PubMed]
- Lin H, Turnbull BW, McCulloch CE, Slate EH. Latent class models for joint analysis of longitudinal biomarker and event process data: Application to longitudinal prostate-specific antigen readings and prostate cancer. Journal of the American Statistical Association. 2002;97(457):53–65.
- Millikan RE, Wen S, Pagliaro LC, Brown MA, Moomey B, Do K, Logothetis CJ. Phase iii trial of androgen ablation with or without three cycles of systemic chemotherapy for advanced prostate cancer. Journal of Clinical Oncology. 2008;26(36):5936–5942. [PubMed]
- Muliere P, Walker S. A Bayesian non-parametric approach to survival analysis using Polya trees. Scandinavian Journal of Statistics. 1997;24(3):331–340.
- Neath AA. Polya tree distributions for statistical modeling of censored data. Journal of Applied Mathematics and Decision Sciences. 2003;7(3):175–186.
- Paddock SM, Ruggeri F, Lavine M, West M. Randomized Polya tree models for nonparametric Bayesian inference. Statistica Sinica. 2003;13(2):443–460.
- Pauler DK, Finkelstein DM. Predicting time to prostate cancer recurrence based on joint models for non-linear longitudinal biomarkers and event time outcomes. Statistics in Medicine. 2002;21(24):3897–3911. [PubMed]
- Pawitan Y, Self S. Modeling disease marker processes in AIDS. Journal of the American Statistical Association. 1993;88:719–726.
- Therneau TM, Grambsch PM, Fleming TR. Martingale-based residuals for survival models. Biometrika. 1990;77:147–160.
- Tsiatis AA, Davidian M. Joint modeling of longitudinal and time-to-event data: An overview. Statistica Sinica. 2004;14(3):809–834.
- Tsiatis AA, De Gruttola V, Wulfsohn MS. Modeling the relationship of survival to longitudinal data measured with error. Applications to survival and CD4 counts in patients with AIDS. Journal of the American Statistical Association. 1995;90:27–37.
- Walker SG, Mallick BK. Hierarchical generalized linear models and frailty models with Bayesian nonparametric mixing. Journal of the Royal Statistical Society, Series B: Methodological. 1997;59:845–860.
- Walker S, Mallick BK. A Bayesian semiparametric accelerated failure time model. Biometrics. 1999;55:477–483. [PubMed]
- Wulfsohn MS, Tsiatis AA. A joint model for survival and longitudinal data measured with error. Biometrics. 1997;53:330–339. [PubMed]
- Xu J, Zeger SL. Joint analysis of longitudinal data comprising repeated measures and times to events. Journal of the Royal Statistical Society, Series C: Applied Statistics. 2001;50(3):375–387.
- Yu M, Law NJ, Taylor JMG, Sandler HM. Joint longitudinal-survival-cure models and their application to prostate cancer. Statistica Sinica. 2004;14(3):835–862.
- Yu M, Taylor JMG, Sandler HM. Individual prediction in prostate cancer studies using a joint longitudinal survival-cure model. Journal of the American Statistical Association. 2008;103(481):178–187.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |