Estimating the Optimal Time to Perform a Positron Emission Tomography With Prostate‐Specific Membrane Antigen in Prostatectomized Patients, Based on Data From Clinical Practice

TL;DR

The paper proposes a Bayesian model to determine the best time for a costly prostate cancer PET scan based on PSA levels and patient history.

Contribution

A hierarchical Bayesian model is introduced to jointly estimate PSA growth and PET-PSMA positivity probability for optimal exam timing.

Findings

The model integrates PSA measurements and PET-PSMA results to estimate optimal exam timing.

This approach improves on current practices by using patient-specific data for better decision-making.

Abstract

Prostatectomized patients are at risk of resurgence, and for this reason, during a follow‐up period, they are monitored for prostate‐specific antigen (PSA) growth, an indicator of tumor progression. The presence of tumors can be evaluated with an expensive exam, called positron emission tomography with prostate‐specific membrane antigen (PET‐PSMA). To justify the high cost of the PET‐PSMA and, at the same time, to contain the risk for the patient, this exam should be recommended only when the evidence of tumor progression is strong. With the aim of estimating the optimal time to recommend the exam based on the patient's history and collected data, we build a hierarchical Bayesian model that describes, jointly, the PSA growth curve and the probability of a positive PET‐PSMA. With our proposal, we process all past and present information about the patients PSA measurement and PET‐PSMA results, in order to give an informed estimate of the optimal time, improving current practice.

Figures

Joint modeling of xi(t) and πi(t). The plot shows the relation between time t and PSA evolution xi(t) (on the left) and between PSA xi(t) and the probability πi(t) of a positive test (rotated, on the right). After choosing the desired probability, the associated PSA and time can be recovered through the model, following the black arrows.

PSA measurements collected on three patients. Patients 1 and 2 are of type BCP but with quite different PSA levels, while 3 is of type BCR.

PSA growth curves (on top) and probability curves (on bottom) for the three patients introduced in Figure 2. Patients 1 and 2 are BCP but with quite different PSA levels, while 3 is BCR.

Individual parameters posterior distributions for the three patients introduced in Figure 2.

Real data application—ROC curves for the results of the simple logistic model and joint logistic model. The AUC for the simple logistic ROC is 0.79, the AUC for the mean joint model ROC is 0.86, while the posterior 95% distribution of the AUC index ranges between 0.78 and 0.86. Results are obtained with R package pROC (Robin et al.).

Real data application—Predictive probability of positive PET‐PSMA results for observed negative (0) and positive (1) examinations of the whole data set (left panel), of misclassified results only (middle panel), and correctly classified results (right panel), obtained through cross‐validation.

Real data application—Posterior mean of the log‐psa (on x‐axis) compared with the observed log‐psa (y‐axis): the black points are estimated through our joint model, and the blue points are estimated through the exponential model. The red line is the bisector. Results are obtained through cross‐validation.

Tables

TABLE 1. S1—Individual parameters. The first column of the table shows the percentage of correctly estimated parameters across the 80 individuals. The second column contains the empirical quantiles at levels [2.5%,97.5%] and the sample mean (the value in the center) of the distribution of R^ across the 80 individuals.

Parameter	Percentage	CIR^
log(λi)	95%	[1.00 1.00 1.01]
τi	95%	[1.00 1.00 1.00]
log(γi)	91%	[1.00 1.00 1.00]
log(μi)	95%	[1.00 1.00 1.01]
log(ai)	100%	[1.00 1.00 1.00]
σi2	96%	[1.00 1.00 1.00]

TABLE 2. S2—Individual parameters. For each scenario (columns) and parameter (rows), the first row reports the percentage of correctly estimated parameters. The second row shows the empirical quantiles at levels [2.5%,97.5%] and the sample mean (the value in the center) of the distribution of the confidence interval (CI) widths. The third row displays the empirical quantiles at levels [2.5%,97.5%] and the sample mean (the value in the center) of R^. All components of the table are computed across the 80 individuals and 100 data sets.

Parameter		Scenario 1	Scenario 2	Scenario 3	Scenario 4
log(λi)	%	95%	93%	94%	95%
	CIW	[3.11 5.51 9.98]	[4.39 7.44 12.76]	[4.11 8.96 18.66]	[3.86 6.73 12.34]
	CIR^	[1.00 1.02 1.10]	[1.00 1.03 1.19]	[1.00 1.04 1.30]	[1.00 1.03 1.18]
τi	%	98%	96%	88%	97%
	CIW	[0.26 1.16 3.52]	[0.40 1.56 4.58]	[0.34 2.25 6.98]	[0.33 1.57 5.07]
	CIR^	[1.00 1.01 1.06]	[1.00 1.02 1.11]	[1.00 1.01 1.10]	[1.00 1.02 1.11]
log(γi)	%	98%	90%	92%	96%
	CIW	[0.10 0.26 0.44]	[0.00 0.28 0.52]	[0.01 0.43 1.93]	[0.08 0.29 0.54]
	CIR^	[1.00 1.02 1.14]	[1.00 1.05 1.33]	[1.00 1.01 1.07]	[1.00 1.03 1.16]
log(μi)	%	95%	90%	88%	93%
	CIW	[0.30 0.72 1.45]	[0.33 0.85 1.67]	[0.4 1.43 3.91]	[0.39 0.80 1.50]
	CIR^	[1.00 1.01 1.08]	[1.00 1.03 1.16]	[1.00 1.03 1.24]	[1.00 1.03 1.16]
log(ai)	%	94%	94%	97%	93%
	CIW	[2.39 4.04 6.36]	[2.72 4.66 8.03]	[2.56 4.52 7.10]	[2.89 4.90 9.12]
	CIR^	[1.00 1.04 1.22]	[1.00 1.08 1.52]	[1.00 1.02 1.10]	[1.00 1.04 1.23]
σi2	%	88%	91%	94%	83%
	CIW	[0.94 1.62 2.89]	[1.39 2.23 3.96]	[1.07 1.95 4.04]	[1.19 2.01 3.78]
	CIR^	[1.00 1.01 1.08]	[1.00 1.02 1.11]	[1.00 1.00 1.03]	[1.00 1.02 1.14]

TABLE 3. Real data application—Global parameters describing the PSA level of the proposed model before the change point. The table shows the limits of the 95% CIs (2.5% and 97.5%), the posterior mean (Mean), and the value of the R^ statistic. Positive coefficients are associated with a lower level of PSA.

Parameter	2.5%	Mean	97.5%	R^
αμ[1] ‐ Intercept	−10.37	−7.76	−5.81	1.03
αμ[2] ‐ Ormono adj	−0.20	1.84	3.84	1.01
αμ[3] ‐ Ormono salvage	−5.67	−3.18	−0.49	1.00
αμ[4] ‐ Radio adj	−1.66	0.58	2.85	1.03
αμ[5] ‐ Radio salvage	−0.53	1.46	3.36	1.07
αμ[6] ‐ Lymphadenectomy	−1.30	1.13	4.25	1.01

TABLE 4. Real data application—Global parameters describing the PSA level of the proposed model after the change point. The table shows the limits of the 95% CIs (2.5% and 97.5%), the posterior mean (Mean), and the value of the R^ statistic. Positive coefficients are associated with a lower level of PSA.

Parameter	2.5%	Mean	97.5%	R^
αγ[1] ‐ intercept	−2.99	−2.58	−2.18	1.01
αγ[2]-PR	−0.14	0.20	0.56	1.18
αγ[3]-PT	−0.29	0.08	0.44	1.05
αγ[4]-PN	−0.78	−0.17	0.47	1.18
αγ[5]-S	0.32	0.64	1.03	1.01
αγ[6] ‐ Age	−0.25	−0.04	0.14	1.01
αγ[7] ‐ Ormono adj	0.18	0.71	1.28	1.13
αγ[8] ‐ Ormono slv	0.26	0.73	1.28	1.02
αγ[9] ‐ Radio adj	−1.69	−1.13	−0.54	1.04
αγ[10] ‐ Radio slv	−0.65	−0.24	0.21	1.03
αγ[11] ‐ Lymphadenectomy	−0.41	0.03	0.49	1.05

TABLE 5. Real data application—Global parameters describing the logistic transformation of the proposed model. The table shows the limits of the 95% CIs (2.5% and 97.5%), the posterior mean (Mean), and the value of the R^ statistic. Positive coefficients are associated with higher probabilities of positive PET‐PSMA examinations.

Parameter	2.5%	Mean	97.5%	R^
αβ[1] ‐ Intercept	−0.39	0.84	2.15	1.00
αβ[2]-PR	−0.10	0.68	1.51	1.00
αβ[3]-PT	−1.04	−0.22	0.62	1.00
αβ[4]-PN	0.06	0.88	1.72	1.00
αβ[5]-S	0.49	1.28	2.09	1.00
αβ[6] ‐ Age	−0.52	−0.14	0.24	1.00
β1	1.64	2.56	3.68	1.00
β2	−0.01	0.00	0.01	1.00

TABLE 6. Real data application—Random effects parameters of the proposed model. The table shows the limits of the 95% CIs (2.5% and 97.5%), the posterior mean (Mean), and the value of the R^ statistic.

Parameter	2.5%	Mean	97.5%	R^
ωμ	2.53	3.14	3.91	1.06
ωγ	0.63	0.75	0.90	1.03
Ψa	−0.41	−0.26	−0.13	1.07
ωa	0.57	0.66	0.76	1.01
aσ2	2.00	2.00	2.02	1.00
bσ2	0.19	0.22	0.27	1.01

TABLE 7. Real data application—Coefficients of the logistic regression used to determine the relationship between administrations of therapies and baseline covariates. Coefficients in bold are associated to p‐values smaller than 0.05. The following abbreviations are used: OA for hormone adjuvant, OS for hormone salvage, RA for radio adjuvant, RS for radio salvage, L for lymphadenectomy, PSA‐s for PSA level at surgery.

Therapy	Intercept	PR	PT	PN	S	L	PSA‐s
OA	−3.22	−0.02	−0.46	1.92	1.20	0.46	1.54
OS	−1.57	−0.08	−0.95	−0.51	0.90	0.54	−1.39
RA	−2.94	2.21	−0.57	0.92	0.35	−0.30	−0.25
RS	−3.14	0.52	1.02	−0.70	1.01	0.23	−0.79
L	−0.03	−0.26	0.84	N.A.	1.12	−0.60	0.25

TABLE 8. Real data application—The table describes the results of the logistic model based on the proposal of Luiting et al. (2020). For each parameter (rows), the first column reports the estimate obtained by maximizing the likelihood. The second column contains the standard error, and the third column shows the p‐value.

	Estimate	Std. error	p‐value
(Intercept)	−0.2754	0.4757	0.5627
PR	0.6082	0.3622	0.0931
PT	0.0918	0.3613	0.7994
PN	0.3738	0.5504	0.4971
S	1.0768	0.3628	0.0030
Age	−0.1381	0.1649	0.4023
log(PSA)	1.3855	0.2628	0.0000
time	0.0070	0.0042	0.0915

Full text

Introduction

1

Nowadays, one of the most important areas of medical statistical development is oncology. The statistical and clinical interest varies from the early detection of tumor presence and locations, the estimation of (personalized) treatment efficacy, the definition of optimal personalized treatment strategies, the analyses of resurgences, the development of ethical trial structures, to many other problems.

Prostate cancer is the most frequent neoplasm in men, with an incidence of around 7% among all new cancer cases International Agency for Research on Cancer IARC, and has therefore attracted a lot of interest in the last 20 years. Research focuses on several topics, starting from the causes and the incidence of the tumor, modeling its growth, analyzing the effect of therapies, and, finally, analyzing the risk, the locations, and time of biochemical relapse (BCR) and clinical recurrences. High levels of serum prostate‐specific antigen (PSA) after primary treatments, such as surgery or radiotherapy, were identified in the literature to be significant indicators of tumor progression taking place at some locations different from the original one (Afshar‐Oromieh et al. 2017; Eiber et al. 2016; Fendler et al. 2019; Hoffmann et al. 2019; Verburg et al. 2016). This means that prostatectomized patients are at risk of developing BCR and metastasis, which is the reason why, during a follow‐up period of several years, they are usually monitored for PCA recurrence, traditionally identified by PSA >0.2 ng/mL.

Positron emission tomography with prostate‐specific membrane antigen (PET‐PSMA) is a nuclear medicine survey, that is currently one of the most sensitive tests for the early detection of tumor presence and location. It is very expensive and complex and, for this reason, patients should be referred to a PET‐PSMA exam only in case of confirmed high risk of BCR. The timing of the exam is very important since if it is performed too late, the patient is subject to excessive risks, while, if it is too soon, there is a high probability of false negative results. The estimation of the optimal time to perform such an examination is still an open problem widely discussed in the literature.

Currently, the most used indicator of BCR is the PSA Doubling Time (PSA‐DT), which is usually monitored over time to control the PSA evolution (Regula et al. 2020). This indicator is calculated using only the last two measurements available for the patient analyzed. According to literature (Regula et al. 2020), a PSA-DT<6 months, together with a thorough clinical examination of the patient stage and conditions, should be used as an indication to perform a PET‐PSMA exam. Some authors model the PSA evolution over time trying to link it to tumor progression (Hirata et al. 2012; Slate and Cronin 1997), distinguishing between patients with and without disease (Carter et al. 1992; Proust‐Lima and Taylor 2009), or accounting of the impact of aging on its evolution (Pearson et al. 1991). Other work mainly focuses on the PET examination (Fossati et al. 2019) or on the correlation between PSA levels and PET results (Pereira Mestre et al. 2019). Finally, some researchers link PSA kinetics and tumor recurrence (Proust‐Lima and Taylor 2009), whereas other authors try to determine the optimal time of PET‐PSMA (Luiting et al. 2020) regardless of the kinetics of PSA. Our work combines many aspects of all these papers into an overall joint statistical model (Duchateau and Janssen 2008; Tsiatis and Davidian 2004; Verbeke and Molenberghs 1997), and exploits its structure to make predictions. While other approaches propose to model either the PSA and the potential resurgences—see, for example, Hirata et al. (2012), Pearson et al. (1991), and Slate and Cronin (1997)—or predict the probability of PET‐PSMA positive results—Fossati et al. (2019) and Luiting et al. (2020)—we propose a joint model that combines both aspects to improve the accuracy of the estimates. We base our model construction on the joint‐model idea proposed by Proust‐Lima and Taylor (2009), in which the authors model the risk of resurgence for a cohort of patients who underwent radiotherapy. They assume that patients after radiotherapy can be divided into classes and that PSA and the risk of recurrence are independent given the latent class membership. The probability that a subject belongs to a class is related to the covariates through a multinomial logistic regression model, with parameters estimated using the maximum likelihood approach. Differently from Proust‐Lima and Taylor (2009), we link the PSA and the results of the PET‐PSMA through a latent description of the true PSA level, which is defined as a patient‐specific functional form. Our aims are also different since, by analyzing a data set from an observational study on prostatectomized patients, we want to develop a methodology that can be generalized to improve the standard of care and be used to predict and provide a measure of uncertainty regarding the optimal time to perform the PET‐PSMA exam.

We assume that the true, yet unobservable, PSA levels depend on patient covariates and on the presence of BCR. Specifically, as in Proust‐Lima and Taylor (2009), the latent PSA level is expected to decrease before BCR and to increase afterward. Both the measured PSA level and PET‐PSMA results are dependent on this latent structure, which is probabilistically linked with them. One of the significant advantages of our approach is that we are utilizing all available information jointly to determine the presence of BCR and the likelihood of a positive PET‐PSMA. Moreover, in contrast to other methods that link PET‐PSMA results to the measured PSA levels, our approach allows for estimating the probability of a positive PET‐PSMA without an associated PSA measurement. By selecting a target probability of having a positive PET‐PSMA and a reliability level, we can estimate the optimal time to perform the test with Markov chain Monte Carlo (MCMC) methods. Inference is performed according to the principles of Bayesian statistics, which allows us to have a complete evaluation of the uncertainty associated with each model component. To do so, a Bayesian network (Chen and Pollino 2012) is built, or more precisely a Bayesian hierarchical model. Bayesian networks are probabilistic graphical models that represent a set of variables and their probabilistic interdependencies, allowing for the data analyst to have an accurate uncertainty quantification of all model unknown parameters. The data set used in this paper is part of a collaboration with the clinicians of the Department of Urology of San Luigi Gonzaga Hospital in Torino (Italy).

The paper is divided as follows: In Section 2, we present the joint model for PSA measurements and PET‐PSMA results. Section 3 is devoted to the definition and the estimation of the optimal time to perform a test. We test our model on simulated data, as explained in Section 4, before applying it to real data: Section 5 contains the results of the model estimated on a group of patients from San Luigi Gonzaga Hospital, in Torino, Italy. The paper ends with some conclusive remarks in Section 6.

Joint Model of PSA Growth Curve and PET‐PSMA Results

2

This section is devoted to explain the statistical joint model: We first describe the growth model for the PSA curve (Subsection 2.1), and then how this model is linked to the logistic structure for the result of PET‐PSMA (Subsection 2.2). Finally, Subsection 2.3 deals with random effects.

The Data and Its Likelihood

2.1

Let xi(t) be the PSA level of the ith patient at time t, where t=0 indicates the time the patient has had a prostatectomy (only operated patients are considered in this work) and i=1,2,…,n, where n is the total number of patients in the study. The recorded variable yi(t) is a noisy version of the non‐negative quantity xi(t), therefore we assume

$$[eqn]$$

The choice of the functional form of xi(t), that is, how PSA levels depend on time and other covariates, is one of the two major components of our proposed joint model (described in Subsection 2.2), but before we discuss it, let us introduce the second component, which is the probability πi(t) for patient i at time t that, if a PET‐PSMA is taken, the result will be positive. For each patient, the test outcome zi(t)∈{0,1} is assumed to be Bernoulli random variable with a patient‐specific and time‐dependent probability, that is,

$$[eqn]$$

Variables xi(t) and πi(t) are not observed, and we can only observe yi(t) and zi(t) at specific time points. The time points where these two variables are recorded can differ within and between patients, and we indicate as Ty,i and Tz,i the set of time points of patient i where yi(t) and zi(t) are recorded, respectively. The set Ti=Ty,i∪Tz,i is the set of points where we have a measure of at least one of the two variables.

Let yio={yi(t)}t∈Ty,i, xio={xi(t)}t∈Ty,i be the observed and latent PSA values over the time points in Ty,i and let zio={zi(t)}t∈Tz,i, πio={πi(t)}t∈Tz,i be the test results and probabilities over Tz,i (i.e., multiple tests for the same patient can be collected, at multiple times, that do not necessarily correspond to times PSA measurements are collected). Conversely, using the superscript u for “unobserved” (as opposed to the previous o for “observe”), let yiu={yi(t)}t∈Tz,i, xiu={xi(t)}t∈Tz,i, ziu={zi(t)}t∈Ty,i, πiu={πi(t)}t∈Ty,i the vectors of variables at the time points where the associated process is not measured, that is, (yi(t),xi(t)) at time points Tz,i where the exams are taken, and (zi(t),πi(t)) at time points Ty,i where the PSA is measured, respectively.

High values of PSA are indicators of tumor progression and are thus associated with a high probability of positive PET‐PSMA results. We will then assume that πi(t) is a function of xi(t), and use this relation to find the optimal PET‐PSMA time since the central idea of this work is to exploit the joint model structure (Proust‐Lima and Taylor 2009; Tsiatis and Davidian 2004) to describe the available longitudinal data. To specify the model, we have to define the joint density of variables (logyi(t),zi(t))′ over Ti, for both observed and missing data. Hence, the missing data are considered further parameters to be estimated during the model fitting, which are easy to estimate within the Bayesian framework.

The joint density for a random sample is then factorized in the following way:

$$[eqn]$$

where f(·) stands for a generic probability density function (to be identified by its arguments) and θ is a vector of parameters. In Equation (3), we are assuming that (logyi(t),zi(t)) are conditionally independent given the latent variables since the connection between the two measurements yi(t) and zi(t) is through the latent PSA level xi(t), and once we know xi(t), no further information is needed to explain the PET‐PSMA result. In the next paragraphs, we illustrate our proposals for f(logyi(t)|logxi(t),θ) and f(zi(t)|logxi(t),πi(t),θ). We want to remark that in Equation (3) we are assuming conditional independence between the measurements, but, on the other hand, as we show in Subsection 2.3, we introduce random effects over some components of θ, to model dependencies across patients.

Joint Model for Each Patient

2.2

To model the time evolution of PSA levels, we assume that it is composed of two phases, the first one, right after prostatectomy, where the PSA level is stationary or even decreasing over time, and the second one, after a patient‐dependent change point in time τi, where it is assumed that the PSA increases, until reaching a plateau indicated as ai. For each patient, the time τi can be interpreted as the unknown time at which resurgence starts, which is a crucial object of inference in our model. The component logxi(t) is modeled as a linear function with patient‐specific intercept λi and slope −μi, for non‐negative μi, that is, if t≤τi

$$[eqn]$$

while for t>τi, we model logxi(ti) using a log‐Gompertz growth model, which is a standard functional form to describe cancer and cancer‐related biomarker evolution (Tjørve and Tjørve 2017; Xu 1987). By the log‐Gompertz model, the evolution function is a weighted mean of the value assumed by logxi(t) at the change point, that is, logxi(τi), and its asymptotic value ai:

$$[eqn]$$

The weight function ei(t) is defined as

$$[eqn]$$

Hence, for t>τi, xi(t) is a version of the log‐Gompertz growth function (Tjørve and Tjørve 2017; Xu 1987) with rate γi∈R+ and asymptote ai∈R+. We assume non‐negative μi and γi to model the decreasing phase before τi and the increasing phase after it.

The second part of the joint model is instead focused on the binary PET‐PSMA results zi(t), and in particular on modeling its probability πi(t). To connect the probability πi(t) to the PSA levels, we use a logistic regression model

$$[eqn]$$

where the linear prediction has a patient‐specific intercept and an extra term linear on time to model temporal trends that are not explained by logxi(t). It should be noted that we define the relation in terms of the true latent PSA levels xi(t), not the observed ones yi(t). We expect β1 to be positive since the larger the PSA, the larger the probability of a positive test. It is easy to see, from Equation (7), that, if β1>0, then πi(t) goes to zero as xi(t) goes to zero:

$$[eqn]$$

All the unknown quantities in the model {(λi,μi,ai,γi,τi,β0,i,β1,β2,σi2)}i=1,⋯,n make up the parameter vector θ.

Random Effects

2.3

We now extend the base model, where each patient has its own set of parameters, to a random effects model. Generally, random effects can help to account for variability and heterogeneity in the data, due to unobserved or unmeasured factors, and may lead to more accurate and reliable estimates of treatment effects.

More precisely, we assume that the model describing PSA evolution can be enriched by the following second‐level distributions:

$$[eqn]$$

where IG indicates the inverse gamma distribution with scale parameter ai,σ2 and shape parameter bi,σ2, and ψi,χ and wχ2 are, respectively, mean and variance of the normal distributions, for χ∈{μ,γ,a,λ}. The log transformations are used to ensure that the parameters are defined in the correct domain. In Equation (9), the means of the normal distributions are allowed to vary from patient to patient since there may exist covariates that affect them via a linear regression as follows,

$$[eqn]$$

where Ci,χ is a patient‐specific vector of covariates of dimension pχ and αχ is a vector of regressors. Finally, covariate information can also be added to the component of the model that is used to define πi(t) by adding, in Equation (7), the following random intercept,

$$[eqn]$$

The random effects (ai,σ2,bi,σ2,ψi,χ,wi,χ2), for χ∈{μ,γ,a,λ}, as well as the hyper‐parameters αχ, χ∈{μ,γ,a,λ,β} are all appended to the parameter vector θ.

Estimating the Optimal Time

3

Equation (7) links the probability of a positive PET‐PSMA test and the latent PSA level; this relation between the two is what we use to find the optimal time. We recall that the PSA level (when observed) is not the true underlying level, but a noisy version of the true latent one that cannot be directly measured.

Defining and Identifying the Optimal Time

3.1

A solution to find the optimal time could be to plug in point estimates of all model parameters (maximum likelihood estimates, or Bayesian posterior means), and then invert Equation (7) to target the desired probability. This strategy is viable, but it does not take into account all sources of uncertainty. The Bayesian approach, which we follow globally to fit the model, gives a better way to estimate the optimal time and, simultaneously, account for uncertainty quantification in a controlled way.

In the Bayesian approach, the overarching goal is to compute the posterior density

$$[eqn]$$

based on which all quantities of interest may be computed. Marginalizing the posterior density in Equation (12), one can obtain, for each fixed time point t, the posterior predictive density

$$[eqn]$$

It should be noted that this can be done since each specific τi is a component of the high‐dimensional parameter vector θ and, similarly, for each t, πi(t) is a parametric function. For each t, we can then verify whether the following condition is satisfied:

$$[eqn]$$

where P stands for posterior predictive probability based on the density in Equation (13), π∗ is a target probability of positive PET‐PSMA, and ρ is a posterior assurance probability (say 95%), similar to a confidence coefficient. It should be noted that in Equation (14) we require t to be greater than the change point τi, whereas π∗ and ρ are design parameters. The resulting decisions are strongly dependent on the latter quantities, which should be carefully chosen in advance, in accordance with clinicians. Default value is ρ=0.95, while π∗ is usually selected with additional receiver operating characteristic (ROC) analysis or according to clinicians expertise. Finally, to select the optimal time ti∗ for the ith patient, we can choose the first available time in the future satisfying Equation (14), that is, the smallest t greater than the largest time in Ti satisfying Equation (14). Figure 1 contains a graphical depiction of the procedure used to obtain the optimal time and how the latent PSA level is connected to the probability of a positive test. In particular, the plot shows the relation between time and PSA evolution, on the left side, and between PSA and the probability of a positive test, on the right side. After choosing the desired probability (on the right‐side x‐axis), the associated PSA and time can be recovered through the model, following the black arrows, from left to right.

Joint modeling of xi(t) and πi(t). The plot shows the relation between time t and PSA evolution xi(t) (on the left) and between PSA xi(t) and the probability πi(t) of a positive test (rotated, on the right). After choosing the desired probability, the associated PSA and time can be recovered through the model, following the black arrows.

Computing the Optimal Time

3.2

Due to the complexity of the hierarchical joint model, the posterior distribution is defined on high‐dimensional data, which prevents us to compute, in closed form, any of the quantities that we may need for inference, as well as normalization constants. As often done with Bayesian models, we use MCMC (Brooks 1998) algorithms to obtain samples from the density in Equation (12) and Monte Carlo integration (Gamerman and Lopes 2006) approaches to approximate the posterior quantities of interest, such as posterior expectations and posterior probabilities and, specifically, the integral in Equation (14), which is the main focus of inference.

Let logyu,b,zu,b,logxu,b,logxo,b,πu,b,πo,b,θb be the bth posterior samples from the associated parameters, where b=1,2,⋯,B, and B is a large number of MCMC iterations. Using these posterior samples, we can approximate the optimal‐time ti∗, since samples from the predictive distribution in Equation (13) are easily obtainable. From Equations (4), (5), and (7), we can see that πi(t), for all t, is a deterministic function of the parameters, even if t is not in Ti. Therefore, for a given t, and each b sample, we can compute

$$[eqn]$$

and

$$[eqn]$$

As a consequence, the set of samples {πib(t),τib}b=1B are from the predictive distribution of Equation (13). It is easy to show that the integral in Equation (14) can be seen as an expectation if we rewrite it as

$$[eqn]$$

where 1·(·) is the characteristic function of a set. Hence, using samples from the predictive distribution, we can approximate the quantity in Equation (17) using a standard Monte Carlo integration, leading to the estimator

$$[eqn]$$

This means that the integral can be approximated by the proportion of posterior samples that are in the set [πi(t)>π∗,τi<t]. The approximation in Equation (18) must be computed for a fine grid of time points and the smallest t≥max(Ti) that satisfies

$$[eqn]$$

is the desired optimal time.

Simulations

4

To investigate the model performance we conduct a simulation study. Due to the complexity of the model, we first examine a single simulated data set and comment on its results before detailing the full simulation study in Section 4.2. The simulation is not intended to mimic a real database; rather, it serves as a test to stress the algorithm and analyze its strengths and weaknesses under challenging conditions. The results of a more realistic scenario can be found in Appendix A.2. In the captions of the tables reporting the results, we will use S1 to indicate results from the single simulated data set discussed in this section, S2 for the results from all data sets, and S3 for the more realistic setting. The results of the model applied to the motivating example, along with a comparison to a competitive approach, are presented in Section 5.

One Simulated Data Set (S1)

4.1

In this section, with a simulated data set, we want to show that the model, and especially the change point values τi, can be estimated using the MCMC algorithm. We simulate a scenario where m=80 patients are undergoing surgery at time 0 and then followed up for several months. For each patient, we simulate the number of elements of Ty,i from the distribution Ud(5,8), where Ud(·,·) indicates a uniform distribution over the integers between the two arguments (included). The times associated with the measurements are sampled randomly without replacement from the set {1,2,⋯,25}, while the numbers of time points of the PET‐PSMA exams are from the Ud(3,5) and the times are randomly sampled without replacement from the set {26,27,⋯,38}. For each i, we assume min{Tz,i}>max{Ty,i}, meaning that the exams are always performed after the last PSA measurement and we have a small set of measurements for each patient. Let Ty,i={ti,j}j=1ni be the set of ordered points: Given these values, to simulate τi we sample from U(ti,3,ti,(ni−2)), so that τi is in the middle of the temporal window. We simulate data in this way to create a challenging situation, where the data points are very few, and the PET‐PSMA exams {zi(t)}i are all performed much later than the last PSA measurement {yi(t)}i. This is done to highlight better how the method can be used to predict the PET‐PSMA results for future time t, for which we have not observed the PSA level, and when data information is poor. For each patient i, we simulate nine dichotomous variables {Cij}j=19, and a final variable Ci10=⌊C^i10⌋, with C^i10∼i.i.d.N(75,7), used to describe the patient age. We then assume Ci,μ=(1,Ci1,Ci2,Ci3,Ci4,Ci5) and Ci,β=(1,Ci6,Ci7,Ci8,Ci9,Ci10). The remaining parameters are

$$[eqn]$$

Note that most of the parameters chosen for the simulation are randomly selected, and there is no intention to be realistic. A Gaussian prior N(0,100) is used on

$$[eqn]$$

and all regression coefficients, where bσ2/(aσ2−1) and bσ22/(aσ2−1)2(aσ2−2) are, respectively, mean and variance of the parameter σi2, while ωa∼IG(1,1). Given the range of value that we expect in the real data application, and the values that we use to simulate the data, the normal priors can be considered weakly informative, as they put most of their mass on likely values without being too informative, and standard choice for the parameters of the inverse gamma are chosen.

The prior of τi is a patient‐specific mixed‐type distribution: τi can assume value in {ti1,[ti2,tiJi−1],tiJi}, where {ti1,ti2,…,tiJi}≡Ty,i, with a probability mass of 1/3 on ti1 and tiJi, and a Uniform distribution in [ti2,tiJi−1]. The reasoning behind this choice is that a change point in [0,ti2) is not identifiable, since only the observation collected at time ti1 is available to estimate the decreasing phase. The same applies in the time window (tiJi−1,+∞) where we could gain information only from tiJi to estimate the log Gompertz coefficients. Since the change point is one of the most important parameters for the final identification of the optimal time, we prefer to have a parameter that is always identifiable. We run our algorithm for 150000 iterations, with burn‐in equal to 100000 and thinning parameter equal to 10. The reason to discard such a long burn‐in sequence is that the model is fairly complex, and we want to make sure that a reasonable hope for convergence holds. In addition, the model contains hyperparameters over which extremely flat prior distributions have been put. As MCMC diagnostics we use the R^ (see Gelman et al. 2014), an index specifically designed to check mixing and stationarity of the posterior chains, by looking at their inter‐ and intra‐variabilities. The algorithm is a Metropolis‐within‐Gibbs MCMC with the adaptive Metropolis steps proposed in Andrieu and Thoms (2008), and a Pólya‐Gamma update (Polson et al. 2013) for (β1,β2,αβ). The model is implemented in R, and for each data set approximately 7 h are needed to obtain posterior samples.

In Table 1, for each individual parameter, we show the posterior 95% credible interval (CI) and evaluate the percentage of parameters correctly estimated (for which the true values used for simulation are inside the 95% CI). From Table 1, we see that the percentages calculated for more than 80 patients are close to 95%, with the median R^ being very close to the optimal value of 1 (Gelman et al. 2014). In terms of global parameters, see Table A1, all CIs of the global parameters contain the true values used to simulate the data, with good R^ values but the CIs of the logistic parameters tend to be large, particularly for αβ[1]. It should be noted that the true value of β2 coincides with the CI lower limit.

Full Simulation Study (S2)

4.2

The full simulation study comprises four scenarios: scenario 1 and scenario 2 differ only in the noise value (determined by parameters aσ2 and bσ2), while the third scenario differs in most of the parameters. Scenario 4 closely resembles scenario 1 but is composed of fewer data points. Specifically, in scenarios 1–3, the data sets consist of 80 patients, with each patient having between 8 and 15 PSA measurements and between 3 and 5 PET‐PSMA examinations. In scenario 4, the data set is composed of 80 patients, with each patient having between 6 and 10 PSA measurements, and between 1 and 3 PET‐PSMA examinations. The underlying parameters used to simulate the data sets are summarized in Table A2. For each of the four scenarios, we simulate 100 data sets and run the MCMC procedure to obtain samples from the posterior. The results are summarized in Tables 2 and A3 and A4, reporting the percentage of individual and global parameters correctly estimated (the true values are within the posterior 95% interval). Additionally, the tables include the quantiles at levels 2.5% and 97.5%, the mean of the interval widths, and R^, computed across the 100 data sets of each scenario.

Except for the variance of the random effects, the model performance in predicting PSA levels tends to be better than in the logistic component. If the number of collected PET‐PSMA measurements decreases, although the percentage of correctly estimated parameters remains relatively high, the amplitude of the posterior credible interval increases, resulting in reduced precision in estimating the parameters for the logistic model. This highlights how the quality of estimation and the reliability of the results strongly depend on the number of available data points for each patient. In particular, the quality of the estimates for the logistic regression parameters depends heavily on the scenario and how much information that scenario provides about the specific parameters, with percentages going as low as 5% (especially, for scenario 3 which turns out to reflect a very challenging case) and, at times, percentages close to 95% being associated with large CI widths. Moreover, the simulation study highlights that higher values of the noise term lead to a decrease in the model efficiency in estimating the growth coefficients μi and γi, which can be seen by the interval widths in Table 2. However, the overall performance of the model is still satisfactory. These results suggest that, for practical application in a clinical setting, a large data set is necessary to accurately learn the relationship between PSA levels and test outcomes. Once this relationship is established, it can be effectively used to predict outcomes for future patients, even in cases where limited data and information are available.

We further assess the model under a more realistic range of parameters, with results presented in Appendix A.2. In this setting, we simulate 100 data sets using the mean parameter estimates obtained from the real case study (detailed below, in Section 5) while assuming a similar data set composition. Under this approach, all diagnostics used to evaluate the model ability to retrieve the true values used in data simulation show improvements, particularly for those parameters that performed poorly in simulations S1 and S2. These results serve as a sensitivity analysis for the case study, confirming the reliability of the model and algorithm.

Application to Clinical Data

5

The database is built on clinical practice in the San Luigi Hospital in Torino. In prostate cancer follow‐up setting post radical prostatectomy (RP), PSA is the main source of information to base clinical decisions regarding interventions such as PET‐PSMA exam. We have n=187 patients and, for each of them, we have several demographic and clinical variables, as well as the PSA measurements and PET‐PSMA results taken at different times after RP. The number of time points in Ty,i ranges from 4 to 17, while the ones in Tz,i from 1 to 4, and the follow‐up period is between 4 and 280 months. We can distinguish between patients who, after prostatectomy, show relatively low values of PSA after RP, but may be eventually subject to a BCR at change point (BCR patients), and biochemical‐persistence patients (BCP patients), who present persistent benign/malignant residual tissue after surgery (conventionally signaled by PSA > 0.2) and are usually treated with therapy to inhibit cancer growth—for them, PSA levels may even initially decrease until a new increase occurs at the change point. The presence of a mixed population of BCR and BCP patients and the possible administration of different therapies make the inference task a very difficult exercise. To present the possible data paths, Figure 2 shows data collected for some anonymous patients, that we refer to from now on as patients 1, 2, and 3. The shape of the data for these patients is quite different (patients 1 and 2 are BCP, while patient 3 is BCR), but also similar shapes can correspond to quite different PSA values (patient 1 and patient 2). Before surgery, each patient is assigned to a different category according to the status (Si), but some more information is collected at surgery time; in particular, pathological stadiation according to TNM classification, and the clinical status of prostate margins. Regarding tumor status PiT, the clinicians use four different categories: T1 (clinically unapparent tumor), T2 (tumor confined within prostate), T3 (tumor extends through the prostate capsule), T4 (tumor is fixed or invades adjacent structures). For the sake of simplicity, here the four tumor categories have been reduced to two: PiT=0 if T1 and T2 and PiT=1 otherwise. The lymph nodes implication is evaluated through a binary variable (PiN), while metastases are always present and, thus, not relevant. Moreover, the clinical stage of the tumor is usually represented using the Gleason Score (Egevad et al. 2002), and we reduced the original nine categories to two: Si=0 if Gleason Score is less than 6 or if it is equal to 3+4, Si=1 otherwise. Finally, prostate resection margins (PiR) are evaluated after the surgery.

PSA measurements collected on three patients. Patients 1 and 2 are of type BCP but with quite different PSA levels, while 3 is of type BCR.

According to clinical evaluation, each patient with BCR and/or BCP can then be administered four different therapies: adjuvant or salvage androgen deprivation hormone therapy (OAi or OSi), adjuvant or salvage radiotherapy (RAi or RSi). Some patients also underwent regional lymphadenectomy (Li). Adjuvant therapies are administered within 6 months from surgery, while salvage ones are performed after 6 months on from surgery. In addition, for each patient, we know the age Ai. All these data are used as clinical and demographic covariates to gain information useful for our model. A descriptive table of the analyzed cohort of patients can be found in the Appendix (see Table A7).

Results

5.1

In particular, we assume that μi, γi, σi2, and ai are random effects, with ψi,a=ψa, ai,σ2 = aσ2, and bi,σ2 = bσ2 (constant for each patient), while the means of both μi, γi depend on covariates, to gain more information from the available data. However, the observed difference between the reported categories of patients, namely BCR and BCP patients, suggests that the parameters λi do not come from a common population. Instead of constructing a bimodal random effect which is out of our scope, we take each λi to be an individual parameter. We estimated three models with different combinations of covariates Cμ, Cγ, Cπ, all suggested by clinical evidence, and among them, we selected the best model using the Watanabe–Akaike information criterion (WAIC) (Gelman et al. 2014), see Table A8 in the Appendix. The results we describe refer to the best configuration selected, namely Model 2. In Tables 3, 4, 5, and 6, we show the global parameter estimates for the best model, with relative 95% CIs.

Covariates Interpretation

5.1.1

Interpretation of coefficients can be difficult and misleading, particularly when referring to therapies. These are the usual difficulties in interpreting causality with observational and clinical practice databases rather than clinical trials. From a clinical perspective, one could expect therapies to decrease the PSA level before change point τ and to reduce the growth speed after it. However, the analysis is not targeting therapeutic efficacy. Therapies are used here as indicators of the severity of the disease rather than for estimating their effect. For this reason, we report the logistic analyses used to determine the relationship between administrations of therapies and baseline covariates (see Table 7), where it is easy to notice that, as expected, patients with the worst clinical situation at the surgery time have higher probabilities of receiving one or more of the analyzed therapies. On the other hand, baseline covariates have relatively simple interpretations. Following the theory, PSA decreasing level turns out to be almost zero for patients that do not receive therapies (see αμ[1]). Gleason score is a significant factor, as higher levels determine an increase in PSA values. In particular, the standard deviation of decreasing coefficient ωμ and increasing coefficient ωγ are quite different. The high level of the former reflects the heterogeneity of the data set: PSA levels right after surgery and preceding the changing points are different for BCR and BCP. Finally, the probability of a positive PET‐PSMA result at time t strongly depends on the PSA level at t (as was to be expected), on the Gleason score, and on lymph nodes implication. In particular, as the PSA level increases, also the probability increases.

Graphical results

5.1.2

In Figure 3, we show examples of the model outputs. In the first row, the fitted curves for the selected patients (first row) are shown, where the gray solid regions are used to show the 95% CI of τi. In the second row, the approximation on the right‐hand side of Equation (18) is computed using thresholds 0.5, 0.7, and 0.9: As expected, the higher the threshold, the lower the probability. Filled circles and squares are used to indicate the true and negative outcomes of the test, respectively. The posterior distributions of individual parameters λi,τi,γi,μi,ai are reported in Figure 4. The posterior distributions λi show differences that reflect the difference between BCP and BCR patients. In particular, it is easy to see that we obtain different fits for patients 1 and 2 of type BCP, in comparison with BCR patients 3.

PSA growth curves (on top) and probability curves (on bottom) for the three patients introduced in Figure 2. Patients 1 and 2 are BCP but with quite different PSA levels, while 3 is BCR.

Individual parameters posterior distributions for the three patients introduced in Figure 2.

It should be noted that, in this case study, data were collected in an observational setup rather than in a clinical trial. Therefore, the timing of measurements was not prespecified but was influenced by the clinician's decisions and the severity of the illness itself, which is a potential source of bias. We recognize that this is a shortcoming of our analysis. However, modeling the timing of the tests, which would require a point‐process modeling approach, would drastically increase the complexity of the model. We do not believe that there is enough information in the data to estimate the parameters of this additional component accurately, but it would be part of future analyses.

Comparison of the Joint Model to the Logistic Model

5.2

We compare the performance of our joint model to the closest method existing in the literature (Luiting et al. 2020), that can be considered the current standard of care. In particular, a logistic regression model is fit with the R function glm (R Core Team) (following the idea proposed in Luiting et al. 2020), including the same baseline covariates as our joint‐logistic model and the observed PSA level. In this logistic model, which we compare to our joint model, the response is 1 if the PET‐PSMA gives a positive result (i.e., the location of the disease is identified) and 0 if not. We apply to the San Luigi Hospital data set presented in Section 5's both methods. Overall, the results obtained from this logistic analysis (see Table 8) are largely consistent with those from our method. Moreover, the significant covariates are also overall consistent with the reference paper (Luiting et al. 2020), based on a different but similar data set. Specifically, we found that tumor status and log‐PSA level have strong and significant effects on the probability of positive PET‐PSMA results, while resection margin and time have weaker effects. Notably, all of these covariates have positive coefficients, indicating that higher values increase the likelihood of positive PET‐PSMA results. To formally compare the two fitted models, we use them to predict the PET‐PSMA binary outputs, on the data set under analysis, and we evaluate the relative ROC curves, shown in Figure 5, computing their areas under the curve (AUCs). Figure 5 shows the ROC curves for the logistic and the joint model, computed using the auc function from the pROC package (Robin et al.). The AUC for the simple logistic model (red line) is 0.79, while the AUC for the mean joint model (black line) is 0.86, with 95% of the posterior AUC values ranging between 0.78 and 0.86, and a median of 0.82. As a result, we can see that the AUC of the simple logistic model is close to the lower quantile obtained using our approach. This suggests that the use of the latent PSA level, which is one of the main points of our proposal, instead of the measured level used in the competing model (Luiting et al. 2020), tends to yield better results. From the mean ROC curves, we can select the optimal π∗ probability, namely the one corresponding to the higher trade‐off between sensitivity and specificity, to be used for predicting new data; it is selected as the point on the ROC closest to coordinates (0,1), which corresponds to π∗=0.63. Moreover, one big advantage of our joint model is that it enables the prediction of probabilities of positive results for future times, for which the PSA level may still be unknown. For example, this prediction of the PSA level and the related probability of positive success can still impact the treatment strategy adopted in the meantime, thus improving the patient benefit.

Real data application—ROC curves for the results of the simple logistic model and joint logistic model. The AUC for the simple logistic ROC is 0.79, the AUC for the mean joint model ROC is 0.86, while the posterior 95% distribution of the AUC index ranges between 0.78 and 0.86. Results are obtained with R package pROC (Robin et al.).

In addition to the logistic regression presented in Luiting et al. (2020) and analyzed in this section, a common procedure currently used by clinicians is the evaluation of PSA‐DT and PSA‐velocity, based on the last couple of PSA measurements. In the literature, PSA‐DT higher than 6 months is significantly associated with a high probability of positive PET‐PSMA results, but no clear and well‐defined guidance is available, resulting in subjective decisions of time to examination, which depend on the clinician's belief and experiences. For this reason, a comparison of the performances of our proposed joint model and this methodology cannot be performed.

A Cross‐Validation Study

5.3

To evaluate the performance of the proposed algorithm on the real data set, we additionally perform a leave‐one‐out cross‐validation (LOOCV) analysis, according to Gilks et al. (1995). We exploit the cross‐validation predictive density sets

$$[eqn]$$

for the PSA levels, where logyito is the set of observed log‐PSA values except for logyi(t), and

$$[eqn]$$

for the PET‐PSMA examination, where zito is the set of observed log‐PSA values except for zi(t).

For each PSA and PET‐PSMA measurement collected, we sample 1000 realizations from the corresponding cross‐validation predictive density, estimated through 150,000 iterations of the algorithm (burn‐in 100,000, thinning parameter 10). To evaluate the goodness of the prediction, we use the accuracy index. The PSA accuracy, that is, the percentage of predictive 95% CI containing the true values of the corresponding measurement, equals 98.20%. On the other side, PET‐PSMA accuracy, that is, the percentage of PET‐PSMA results correctly predicted (for each data, the threshold is selected through LOOCV on ROC curves), equals 72.24%. Moreover, to validate the posterior distribution of πi(t), Figure 6 shows the probability of positive results across different subsets (0 for negative and 1 for positive) PET‐PSMA results. In the left panel, we present the overall data set, the middle panel focuses on misclassified instances, and, finally, the right panel highlights correctly classified cases. The plot shows that the posterior median of πi(t), respectively, for the outcome 0 and 1, is approximately 0.38 and 0.82 when considering the entire data set. However, when we restrict our analysis to the correctly characterized data, these values shift to ≈0.33 for 0 and ≈0.85 for 1, which demonstrates that the right decision can be made with a high degree of confidence. In the middle panel, when the method encounters challenges in classifying the data, there is a notable increase in uncertainty. The values of πi(t) are close to the threshold π∗ (potentially around 0.55), confirming that the methods fail when uncertainty is high. This underscores the sensitivity of the model to situations with increased ambiguity and suggests that the uncertainty in predictions increases when the method encounters difficulty in making accurate classifications. Moreover, we compare the performance and the accuracy of our joint model with results obtained with standard methods, consistently with existing literature practices. We conduct LOOCV to assess the performance of the PET‐PSMA logistic model—based on Luiting et al. (2020) and presented in Section 5.2—but also the exponential‐fit model, a simpler method that predicts the future (and yet unknown) PSA level at time t3 based on pairs of PSA measurements (PSA1,PSA2) collected at previous times t1 and t2, utilizing an exponential fit (Benecchi et al. 2008; Memorian Sloan Kattering Cancer Center MSK). Note that the PSA level at the given time t3, is predicted by only using two previous recorded values and not the entire history. While this technique offers predictions, its efficacy is limited compared to our joint model, as shown in Figure 7. It illustrates the difference between true and predicted log‐PSA levels, estimated with both exponential and joint models. The plot distinguishes predictions made with the exponential model (1647 predictions shown in blue) and our joint model (1671 predictions represented in black), underlying the limited predictive capacity of the exponential basic model, which handles 24 fewer data predictions due to the intrinsic structure of the prediction mechanism (points referring to time t1 and t2 can not be handled). Nevertheless, it is worth noting that the variability of the predicted values with the exponential model is higher than that with the proposed joint model. This highlights the superior predictive performance of the proposed joint model. Finally, with a standard LOOCV approach, the simple logistic model achieved an accuracy of 65.42% for predicting PET‐PSMA examination outcomes (with a balanced accuracy of 66.27%). In contrast, our joint model improves on these results, reaching an accuracy of 70.83% and a balanced accuracy of 72.24%. We included balanced accuracy as an additional evaluation metric to address potential errors arising from the imbalance between 0 and 1 PET‐PSMA results in the data set, which may be influenced by biological factors; see Stephenson et al. (2006) and Kupelian et al. (1996).

Real data application—Predictive probability of positive PET‐PSMA results for observed negative (0) and positive (1) examinations of the whole data set (left panel), of misclassified results only (middle panel), and correctly classified results (right panel), obtained through cross‐validation.

Real data application—Posterior mean of the log‐psa (on x‐axis) compared with the observed log‐psa (y‐axis): the black points are estimated through our joint model, and the blue points are estimated through the exponential model. The red line is the bisector. Results are obtained through cross‐validation.

Final Remarks and Conclusions

6

Correct and quick identification of the locations of possible metastasis in prostate cancer is a challenging open problem. Despite the availability of several new techniques, their calibration is still debated. In particular, the results obtained using the sensitive nuclear examination known as PET‐PSMA can be improved if correctly combined with a good estimation of the optimal time to perform the examination. The previous sections contain our proposal on how to estimate the optimal time to perform PET‐PSMA, which exploits information from the whole history data of each patient. We have introduced the joint model approach, addressing both PSA growth and the probability of a positive PET‐PSMA, enriched with random effects to enable predictions for future patients. Our proposal is therefore not just a method for predicting individual PSA growth curve and time to PET‐PSMA, but a proposal to drive optimal decisions regarding the patient, which is the core of personalized medicine. In the paper, after explaining the model structure and the joint approach, we have discussed the optimal time estimation. Finally, the model has been estimated both on simulated and real data. Simulations were used to test the proposed model under challenging settings. In particular, we showed that the resurgence changing point time τ is difficult to estimate, as well as the regressive parameters entering the mean of the random effects. On the other hand, simulations also highlight the adaptability of the method to quite different growth patterns of PSA. The results obtained on real data, for an optimal probability π∗ selected via LOOCV, with a confidence level of ρ=95% on the result, give new insights into the model applicability and performance. However, the timing of the PET‐PSMA examinations not being scheduled a priori is a possible source of bias, as it depends on the natural history of the disease. Although ignoring this component is a common practice, we are currently considering a possible extension of our model that can address this issue if the data and covariates contain enough information. On the other hand, both the growth model estimated patterns and the logistic results are easy to interpret and align with clinical evidence. Further research directions could include the implementation of a graphical interface to help clinicians easily exploit the model: Any new patient under analysis should be easily included in the database through the interface, enforcing the model accuracy, and enabling good and quick predictions.

Author Contributions

M.A. preprocessed the data, developed the model, implemented the code, and wrote and revised the paper. G.M. developed the model, implemented the code, and wrote the paper. M.G. formulated the model and revised the paper. S.DL. provided with the data, helped with the clinical interpretation of the results, and revised the paper.

Conflicts of Interest

The authors have declared no conflicts of interest.

Supporting information

Supporting Information