Our data consists of 1,028 women diagnosed with breast cancer in Lahore, Pakistan. All women observed recurrence between February 2011 and February 2018 after initial treatment. They were treated at the same hospital (Institute of Nuclear Medicine & Oncology Lahore, Pakistan). The primary endpoint of this study is death due to breast cancer. Exclusion criteria were: incomplete or missing information, women diagnosed with another disease or another cancer before breast cancer, and bilateral carcinomas. Women who were still alive (survived) at the end of the study, or died due to another reason than breast cancer, are considered right-censored. Age at diagnosis, age at recurrence, estrogen receptor (ER), progesterone receptor (PR), human epidermal growth factor receptor 2 (Her2), tumor grade, radiotherapy and chemotherapy are the predictors included in this study, all were chosen with the help of clinicians and oncologists.

In the data, age at diagnosis and age at recurrence (in years) are continuous variables, estrogen receptor, progesterone receptor and Her2 are represented by binary variables (0 = Negative, 1 = Positive), tumor grade is represented categorically ( I , II , III ). In addition, chemotherapy and radiotherapy are indicated by dummy variables (0 = No, 1 = Yes), here 0 indicated the patients who did not receive the treatment, while 1 = Yes meant they received treatment. Survival time was considered from recurrence to death or drop out, and censoring status was coded with 0 for censored, and 1 for death due to breast cancer. Along with the aim of this study, which is the comparison of different parametric models, it is also a major interest to find out how the two treatments (radiotherapy and chemotherapy), in combination with the other predictors, affect the survival of breast cancer patients after recurrence.

Proportional hazards have been checked by scaled Schoenfeld residuals (Extended data). Furthermore, the statistical tests revealed that not all covariates are statistically significant (p < 0.05), but the global test is statistically significant. Therefore, we can assume that the proportional hazard assumption holds.

In this work, right censored time-to-event data are considered. Survival time is denoted by T and censoring time by C j , where j = 1 , 2 , 3 , … , m are women diagnosed with invasive breast cancer who observed first recurrence after primary treatment. The j th event time is defined by t j = min T j C j , while time for censored and uncensored events is denoted by δ j = I T j ≤ C j . The general relationship form, between T and X is given as T = f X + ε , (1) where X = x 1 , x 2 , x 3 , … , x p , is the vector of predictors which may have an impact on time to failure, and f is the unknown functional form, which may be linear or non-linear. Practically, a 100% exact true relationship is not possible. To cover up uncontrolled chances of error, ε is also included in the model. S t = P T > t is the survival function, which represents the probability of a woman survivor up to a time point, and the hazard rate is written as π t = lim ⍙ t ⟶ 0 P t ≤ T < t + ⍙ t | T ≥ t ⍙ t (2)

The hazard rate represents a probability of an instantaneous failure per unit of time given that an individual patient has survival after time t . ⁵ The main objective of time to failure studies is to estimate the hazard function accurately. For this purpose, different modelling approaches are applied.

In the multivariate approach, the semi-parametric Cox PH model is most popular for the analysis of time to failure data. The general form is written as π t x j = π 0 t exp x j T α , j = 1 , 2 , 3 … , m (3) where, m is the number of patients under study, and x j = x j . 1 x j . 2 x j . 3 … x j . p T is a row vector of p predictors for subject j , while α is the vector of regression coefficients. ⁷ The partial likelihood method employed for estimating unknown parameters is suggested by Cox. ¹⁶ In time to failure analysis, continuous predictors are often categorized, which has disadvantage of information loss within categories. Useful available statistical methods of handling continuous predictors are fractional polynomials (FPs) and restricted cubic splines (RCS).

FPs provide flexible parameterization of continuous predictors, were first used for modelling families of curves by Royson and Altman, ¹⁷ with polynomial of degree m , it can be written as FP m X = α 0 + ∑ j = 1 m α j f j x 1 (4)

m is a positive integer and α 0 , α 1 , α 2 , … , α m are regression parameters, while f j x 1 is divided into two parts: one consists of x 1 p j when p j ≠ 0 , and another of ln x 1 for p j = 0 . ¹⁷ FP models have a wide variety of shapes based on different transformations. The main issue with fractional polynomials is in choosing a suitable power for polynomials, as this has a direct positive relationship with flexibility. To increase flexibility, a greater power can be used, but with the major threat of non-locality. That means that the fitted function at a given point of x 0 depends on data points which are very far from that reference point. ¹⁸

Spline regression is an improved technique, used to overcome the non-locality of fractional polynomials. In spline models, the dataset is divided into multiple parts and these parts are joined with knots. In time to failure regression analysis, spline modelling is extensively used to smooth non-linear effects of continuous predictors. Spline f X has a smooth function. The major problem occurs in choosing the number of knots, as no hard and fast rule is available to apply the suitable number of knots. ¹⁹

Under RCS, the best way of choosing the number of knots is using the quotient of the difference between the largest and average uncensored log survival time and the largest and smallest uncensored log survival time. ²⁰ Royston suggested another method in this respect; according to him, a good way to choose the suitable number of knots is to randomly apply different number of knots every time, and select the best model with measure of information criteria. ²¹ Flexible parametric models have scaling for proportional hazards or proportional odds, which are usually based on transformation of survival function by a link function ή S t x j = ή S 0 t + α T x j (5)

S 0 t is the baseline survival function, α is the vector of unknown parameters for predictors x . Piecewise exponential models (PEMs) are also a reasonable approach to estimate hazard ratios more accurately. ²¹ Under the PEM modelling approach, follow-up time is divided into i intervals, by assuming a constant baseline hazard in each interval, so that π 0 t = π i simplifies to, π j t x j = π i exp x j T α (6)

For the cut points of follow-up times, minimum to maximum time is divided into finite intervals, 0 = n 0 < n 1 < n 2 < n 3 < … < n i = t max . Here, n 0 to n i are time intervals, and the hazard rate of exponential distribution is constant over time intervals. Censoring and all unique event times can be used as time interval cut points, but no hard and fast rule is available to choose cut points, which is a point to ponder; too small or too large cut points may cause under- or overfitting.

One approach to deal with this cut point problem is an extension of PEM, in which a large number of cut points are used and the hazard is estimated semi-parametrically. This is called the Piecewise exponential additive mixed model (PAMM), in which a hazard is modeled through a smooth nonlinear function. In PAMM, predictors contribute to the hazard additively, imposing a quadratic penalty on the basis coefficients: π j t x j = exp f 0 t j + ∑ n = 1 p f n x j . n t j , ∀ t ϵ n i − 1 n i (7) where, f 0 t j denotes log baseline hazard rate, and f n x j . n t j represented effects of smooth nonlinear constant predictors, while t j is finite time cut point. PAMM is an extension of PEM, in which modeling is done by using baseline hazard as a spline basis function, hazard is constant across intervals, through penalization over fitting is avoided even with large number of time intervals. ²²

In fully parametric PH modelling, baseline hazard function is assumed to follow a specific distribution and coefficients are estimated via maximum likelihood. A number of different parametric PH models are derived by applying distributions, such as exponential, Weibull, and Gompertz.

An alternative to parametric PH is AFT models, the corresponding log-linear form of the AFT model with respect to time is given as log T j = α 0 + ∑ j = 1 p x j ′ α j + W ε j (8) where, α j is the vector of coefficients of unknown parameters, W is the scale parameter, and ε is the random error term. ²³ AFT models measure the direct effect of predictors on the survival time rather than the hazard. ε j , as a random variable, assumes different distributions for survival time T , such as enponential, Weibull, gamma, generalized gamma, log-normal, log-logistic and generalized F. ²⁴

The likelihood estimates are maximized using the Newton Raphson procedure, ¹⁵ which may be time consuming and tricky without computer programming. The freely available R software is used to implement all modelling techniques.

Comparison of fitting models is done via measures of fit, which describe accuracy of fitted models for a given data set, usually called goodness of fit measures. Model fitting accuracy has nothing to do with the predictive ability for external data prediction. The Akaike Information Criterion ( AIC) ²⁵ and the Bayesian Information Criterion ( BIC) ²⁶ are two of the most common measures which are used to compare models’ performances. For the PH model AIC and BIC are based on log-partial likelihood AIC = − 2 Log − partial likelihood + 2 df (9) BIC = − 2 Log − partial likelihood + n df (10) where, df represents degrees of freedom of the fit, and n is the total number of observations in the data. The minimum value of both AIC and BIC is considered a good one. Basically, the AIC criterion is used as a penalized function, as if one adds a variable, sampling variability also increases. While BIC imposes stronger penalty in the inclusion of additional covariates to the model. Hurvich et al. ²⁷ suggested a modified AIC, in which n × df + 1 / n − df + 2 is used as a penalty term. A corrected version of BIC proposed by Volinsky and Raftery ²⁸ can also be used, which replaced n with uncensored observations. Corrected AIC ( AICc) and BIC ( BIC c) have written forms as AICc = − 2 Log − partial likelihood + 2 n df + 1 n − df + 2 (11) BICc = − 2 Log − partial likelihood + n uncensored df (12)

According to the Ethical Guidelines for Epidemiologists (IEF-EGE) and the regulations of the ethics committee located at the Advanced Studies and Review Board, University of the Punjab, Lahore (Pakistan), no ethics approval is needed, because the analysis is based on routine data. The study was critically cleared by the Advanced Studies and Review Board of Punjab University. The letter of support written by the departmental head was submitted to the selected hospital. Prior to data collection, written consent was obtained from the head of oncology department and confidentiality was maintained by coding from data collection to analysis.

There are several limitations to our study. First, the use of a single case study may be viewed as a limitation. However, a simulation study can also be designed to validate results. Second, the model comparisons are based on within sample information measures ( AIC, AICc, BIC, BICc), while predictive performances of models can also be checked via different measures. Third, sensitivity analyses of choosing different numbers of knots in spline-based models can be performed to make firm conclusions.

Data is available from the corresponding author, Dr. Florian Fischer ( florian.fischer1@charite.de), upon reasonable request. Data can be used for research purposes, but cannot be published because it is taken from a hospital.