Current study focuses on 1504 men with primary PCa and treated by EBRT, ADT, prostatectomy, or a combination of these interventions. The dataset for this study was obtained from a renowned cancer hospital in Pakistan. Ethical approval was secured from the departmental head and dean of sciences at the University of the Punjab, with permission granted by the hospital authority. Patients registered at the hospital between 2012 and 2019, diagnosed with PCa, and followed for at least two visits were included in the study. The data were manually entered into Excel sheets from patients’ files, ensuring completeness and accuracy while excluding any insufficient information. During follow-up, treatment efficacy is assessed by observing the patients’ tumor status, detecting local or distant recurrence, and noting instances of death, whether attributed to PCa or not. Specifically, our focus is on the time taken for tumor shrinkage, a parameter monitored by physicians through regular visits and observing PSA measurements after the initial treatment. PSA measurements were collected between the end of initial treatment and the occurrence of study event, which is tumor shrinkage up to satisfactory level or event not happened at the end of study period. As shown in Table 1, this study observed post-treatment PSA levels in the log scale (logPSA), BMI (kg/m2), time to tumor shrinkage (in months), and treatment (categorized into 7 groups: 2, 3, 4, 5, 6 versus 1 = ADT as a reference) variables.

Repeated PSA measurements are taken during check-ups; its increased level after treatment indicates the growth of cancer cells. The PSA individual trajectories collected between the end of initial treatment and the occurrence of the event are depicted in Fig. 1. In general, this longitudinal process highlights variations in the biomarker’s long-term changes, reflecting both “tumor shrinkage” and instances of “censorship”, for example, based on study aim which is to illustrate potential relationship between PSA and time until patients achieve a satisfactory tumor status as directed by physicians.

Individual and mean profiles of observed PSA data over time

The time to tumor shrinkage represents the outcome variable categorized as ‘Yes’ for patients with the event of interest or ‘No’ for those patients who did not experience tumor shrinkage or left the follow-up study. ‘Yes’ is coded as 1 and ‘No’ is coded as 0.

Considering the impact of PSA levels on PCa patients’ recovery, utilizing a statistical model is crucial for understanding the relationship between PSA measurements and tumor status. The dynamic progression of PCa varies among patients, highlighting the PSA biomarker’s significance in describing disease progression and its’ correlation with tumor status. This potential is unveiled through the combined analysis of repeated PSA measurements and time to tumor shrinkage variables. Table 1 illustrates the baseline characteristics of PCa patients. The median follow-up number of times is 3 per patient with a range of 1 to 5, which are distributed unequally among individuals. Two outcomes are distributed as \(log PSA \left(1.96\pm 2.03\right),\) and shrinkage of tumor (Yes, No). The event of interest for this study is individuals’ condition (1: tumor shrinkage, 0: right censored) at the end of follow-up time, from 1,504 patients 960 observed events of interest, and 544 were right censored.

A joint model [21, 22], incorporating mixed-effects and event time components, has been developed to capture the relationship between PSA and tumor status. The mixed-effects model describes the evolution of PSA over time, taking into account both fixed and random effects. The event time sub-model, employing the Cox proportional hazards model, analyzes censored data, with baseline predictors including BMI and treatment.

A simple mixed-effects model for longitudinal data with general form is written as,

$$_\left(t\right)=_\left(t\right)+_\left(t\right),$$

(1)

where, \(_\left(t\right)\) is mean of predictor for both fixed and random effects, \(_\left(t\right)\) is an error term. Investigation of the effect of covariates on repeated PSA measurements can be applied using quadratic, cubic, or non-parametric fits.

An event time sub-model is formulated for censored data, as not all patients experienced an event of interest. Mostly right censoring occurs, due to dropouts before the end of follow-up time [23]. Cox proportional hazards (PH) [24] is the most popular semi-parametric model to analyze event time data [25], which is formulated as,

$$_\left(t\right)=_\left(t\right)\text\left(^_\right),$$

(2)

where \(_\left(t\right)\) is a baseline hazard function at time \(t\), \(_\) denotes baseline predictors for regression coefficients’ vector \(\gamma .\) The Kaplan-Meier event curve in Fig. 2 illustrates the time to low-status tumor following treatment combinations. It indicates an increased probability of time to low tumor with the administered treatments, as compared to other treatments ADT observed more effective in terms of time to low tumor status. The baseline hazard can either be unspecified or parametrically modeled using appropriate distributions such as Weibull, Gamma, Exponential, and others [26].

Kaplan-Meier estimates of the probability of survival for individuals on each treatment

The joint model specifies the hazard of the event, which is dependent on individual characteristics of its longitudinal trajectory, as follows

$$_\left(t\right)=_\frac\left\_^+\varDelta t|_^\ge t,}_\left(t\right),_\right\}},$$

(3)

where, \(}_\left(t\right)=\_\left(s\right),0\le s<t\}\) is a history of unobserved longitudinal process \(_\left(s\right)\) up to time t, and \(_\) is a vector of time-varying covariates.

Constructing a joint model involves integrating various association structures to unify longitudinal and event time processes. Commonly used association structures include current value, shared random effects, and current value and slope [27]. The current value association structure assumes that the true value \(_\left(t\right)\) of longitudinal measure at time \(t\) is predictive of the risk of experiencing an event at the same time. The Cox’s PH sub-model with this association structure is written as,

$$_\left(t\right)=_\left(t\right)\text\left(^_+__\left(t\right)\right),$$

(4)

\(\alpha\) is a vector of associated parameters to quantify the association between longitudinal process and hazard for the event at time \(t\). It is interpreted as one unit increase in current value is associated with \(\text\text\text\left(_\right)\) increase in risk of event at the same time, given that event has not occurred before \(t.\)

In a shared random effects association structure, random effects from the longitudinal sub-model are incorporated into the relative risk sub-model as linear predictors, facilitating the sharing of random effects between the two [28]. Another association structure, current value and slope, establishes a linkage between event time and longitudinal sub-models by adding the rate of change of measurement at time \(t\) estimated by taking the derivative of \(_\left(t\right)\) with respect to time. A sensitivity analysis is conducted to select the appropriate association structure. The structure is chosen based on the BIC criteria, opting for the model with the lowest BIC value.

For dynamic prediction, the three-step process is employed to develop a prediction model that accounts for both baseline patient characteristics and longitudinal measurements of PSA values.

We utilized a series of joint logPSA and time-to-tumor-status models, exploring diverse options for both longitudinal and event time sub-models, along with various association structures for joint modelling. The first step describes the evolution of PSA measurements over time, and the second step utilizes this information to model event time. Finally, dynamic prediction is performed using the proposed joint model. As a preliminary step, a covariate selection process is carried out for a sub-model of the longitudinal outcome, and heterogeneity in residual plots is mitigated using a logarithmic scale of PSA.

A mixed-effects model is proposed for the evolution of PSA over time to account for the positive correlation between observed measurements within the same patient. This model includes time (in months) and baseline treatment variables (ADT, prostatectomy, EBRT, and combinations of these). Based on BIC [29] criteria, the best model for logPSA repeated measures is formulated as,

$$_=(_+_)+_+_)_+__^+__+__\times _+__\times _^+_\left(t\right)$$

(5)

In the event time process, following the initial covariate selection, treatment and BMI are identified as significant covariates.

$$_\left(t\right)=_\left(t\right)\text\text\text(__+__),$$

(6)

Joint model prediction based on event time probabilities and future PSA observations for any new patient \(j\), utilizing longitudinal PSA values \(}_\left(t\right)=\_\left(s\right),0\le s<t\}\) and baseline covariates \(_.\) Conditional probability \(_\left(u\right|t)\) is used to predict about patient \(j\), who will have low tumor for time \(u>t\) observing PSA. Using information available at time \(t\), prediction is updated dynamically at any follow-up visit time \(^},\) such that \(^}=t<^}<u\) to produce a new prediction \(_\left(u|^}\right)\).