As the second leading cause of death globally, cancer was responsible for almost 9.6 million deaths in 2018 (1). According to the reports of the World Health Organization (WHO), about 1 in 6 deaths worldwide is due to cancer and approximately 70% of cancer deaths occur in low- and middle-income countries (2). Breast cancer (BC), with almost two million new cases in 2018, is the most common cancer in women and the second common cancer after lung cancer in the world, leading to the death of 626,679 BC patients (1). However, the cure rate of BC will be high if it is quickly diagnosed and treated according to the best available therapies (3). In the Islamic Republic of Iran, BC constitutes more than 24% of all cancers (4) with an incidence rate of 24.8 to 34 per 100,000 women and a mortality rate of less than 10.2 per 100,000 women in 2018 (1). This cancer affects most Iranian women between the ages of 35 and 44 years (5, 6), which is more than 10 years lower than the age of BC onset among women in the Western countries (4). Based on the results of previous studies (7), BC in Iranian women is often diagnosed at later stages; hence, further treatment is barely responding, and the patients have lower survival rates.

Survival analysis investigates prognostic factors of survival in patients using methods such as the Kaplan-Meier method and the Cox proportional hazards model (8). The Cox model is the most common statistical model in analyzing survival data (9). The term “proportional hazards” means that the risk for a patient is proportional to the risk for another patient and it does not change over time (10). In this model, the regression coefficient is estimated by maximizing the partial likelihood (PL) function, and the baseline hazard is an unspecified positive function of time that determines the shape of the survival function. It is not often estimated, and the positivity constraint on it is not considered, leading to the loss of efficiency.

Due to its importance, researchers are still looking for an optimal estimate of the baseline hazard function. In recent years, a new method has been proposed by Ma et al. (11) to simultaneously estimate the baseline hazard function and regression coefficients, which is called the maximum penalized likelihood (MPL) method. Although the MPL methods have already been existed (12-14), they have deficiencies such as not ensuring the positivity constraint on the baseline hazard. In the method of Ma et al. (11), a new iterative optimization algorithm combines Newton’s method with a multiplicative iterative algorithm to estimate the baseline hazard function, thus respecting the positive constraint on it. Moreover, this algorithm provides accurate variance approximations for both regression coefficients and baseline hazard (11, 15, 16).

Because of the unprecedented growth of BC, especially in developing countries, it is necessary to use new efficient methods for analyzing predictive factors of this disease. Therefore, this study was conducted to determine the survival rate and the risk of death in BC patients and to estimate the survival predictive factors using the new MPL method. Moreover, a comparison is made between this new method and the usual PL method.

3.1. Data Source

This prospective cohort study used the data of 2300 patients with BC registered at the Cancer Research Center of Shahid Beheshti University of Medical Sciences in Tehran, Iran. Established in 2006, this governmental, referral, and comprehensive cancer control center has countrywide admission. It has a linked clinical section, Comprehensive Cancer Control Center, which focuses on early detection, diagnosis, and treatment of different cancers.

The survival outcome of 234 patients was not available in the dataset. Thus, they were excluded from the study. Of 2066 remaining patients, 1710 patients were excluded because of incomplete medical records, male gender, and death due to a cause other than breast cancer. Therefore, the final analysis was done on 356 patients diagnosed with BC between 1999 and 2015. This sample size was sufficient based on Kleinbaum and Klein 10, considering a recurrence variable and using the Equation 1:

Equation 1.

Where α = 0.05 and β = 0.20 are probabilities of type I and type II errors, Δ = 1.6 is the effect size, and P_EV1 = 0.75 and P_EV0 = 0.10 are the probabilities of death in the recurrence and non-recurrence groups, respectively. Thus, the power was 0.80, and the sample size was obtained as 347. Two observers examined the medical records of patients with the kappa coefficient of 0.87. Figure 1 represents the flow diagram of data extraction.

Figure 1. The flow diagram of data extraction

3.2. Variables and Statistical Methods

In this study, the survival time of patients (months) and predictive variables were taken into account (Table 1). First, the Kaplan-Meier curves were drawn for all predictive variables to peruse the survival of patients at different levels of predictive variables. The log-rank test was used to conduct univariate analysis and make comparisons based on survival at different levels of predictive variables to find factors with a significant effect on patients’ survival time. By extracting significant factors in univariate analysis, multiple analysis was done using the Cox proportional hazards model. The Cox model is written as Equation 2

Equation 2.

where X_i is the vector of predictive variables for the ith patient, β is the vector of coefficients, and h₀(t) is the baseline hazard that is an unknown and non-negative function of survival time (t) with no assumptions about its shape (17).

Regression coefficients of the proportional hazard model are usually estimated by maximizing the Cox partial likelihood (PL) function (18) where the baseline hazard is not required when estimating the regression coefficients. MPL is another new method that considers the simultaneous estimation of baseline hazard and coefficients in the Cox model. In this method, optimization is achieved using a new iterative algorithm, which combines Newton’s method with the multiplicative iterative algorithm (19) and satisfies the non-negativity requirement on the baseline hazard estimate (11).

The Cox model depends on an assumption that if it is not valid, the analysis will silently give misleading results, in the sense that the size of the effects or even their direction may be inaccurate. For this reason, the validation of the assumption is done by the scaled Schoenfeld residuals (20).

The correlation between variables was also investigated in order to avoid multicollinearity. Multicollinearity is a condition in which some of the independent variables are highly correlated, causing large standard errors of estimated coefficients. In order to reduce the standard errors of the estimated coefficients, one proposal is to drop one or several independent variables from the model (21). The significance level was set at 0.05 and the data were analyzed by the survival MPL package of R version 3.4.3 software.

Of the 356 studied patients, 66 (18.5%) died by the end of the study and 290 (81.5%) patients censored. The age of patients ranged from 24 to 84 years with a mean ± standard deviation (SD) of 48.42 ± 11.27 and a median of 48 years.

Figure 2 shows the Kaplan-Meier survival curve for the variable “stage of cancer” (due to the lack of space, the Kaplan-Meier curve is presented only for this variable). Based on the curve, patients with stage I had higher survival than patients in other stages. Table 1 shows the clinical, pathological, and biological characteristics of the patients. The last column of Table 1 presents the results of the log-rank test. According to the P values obtained from this test, there was a statistically significant difference between the levels of Kaplan-Meier survival curves in terms of variables including the stage of cancer, tumor size, tumor grade, lymphovascular invasion, cancer recurrence, estrogen receptor, and progesterone receptor. Therefore, only were these variables used in multiple analysis.

Figure 2. The Kaplan-Meier survival curve for cancer stage

Multicollinearity between the variables and proportional hazard assumption were first checked. The highest correlation (0.78) was established between the stage of cancer and tumor size. Thus, the stage of cancer was not included in the multiple analysis. In assessing the assumption of proportional hazards, the P values obtained from the correlation check between Schoenfeld residuals and the ranked failure times for any of the variables were not significant, indicating that there was not enough evidence to reject the assumption of proportional hazards. Table 2 shows the results of the Cox model fit using the PL and MPL methods.

Based on the results of Table 2, for most variables, the standard errors were lower in the new MPL method than in the PL method, and the estimation of the parameters and the hazard ratio were different from those of the PL method. For example, the standard error and hazard ratio were 0.727 and 47.87, respectively, for the recurrence variable in the PL method while these values were 0.621 and 33.21 in the new MPL method.

In multiple analysis using the PL method, statistically significant relationships were found between survival and lymphovascular invasion (P = 0.01) and recurrence (P < 0.001). Based on the new MPL method, in addition to lymphovascular invasion (P = 0.007) and recurrence (P < 0.001) variables, estrogen receptor (P = 0.045) was also a significant variable. Based on these results, the risk of death was 1.72 times higher in patients with negative receptor variable than in patients with a positive receptor variable (HR = 1.72).

Figure 3 shows the plot of the overall survival for patients with BC with a 95% confidence interval based on the new MPL method and the estimated baseline hazard.

Figure 3. The plot of overall survival for patients with breast cancer

Table 3 shows the patients’ survival rates of 12, 24, 60, 84, and 120 months using the nonparametric Kaplan-Meier method and the MPL method by estimating baseline hazard function. According to Table 3, the new MPL method estimated the survival probability after 24, 60, 84, and 120 months with fewer standard errors compared to the nonparametric Kaplan-Meier method.

Cancer is the second leading cause of death globally and it was responsible for almost 9.6 million deaths in 2018 (1). BC is a malignant tumor mostly observed in women. BC affects Iranian women at least one decade sooner than women in developed countries (22). Reports show that the average age in BC occurrence is 61 years among US white females (23) while our study, like other studies from Iran (24, 25), showed that the mean age at BC diagnosis is about 48 years. To estimate survival in patients with breast cancer, researchers mainly employ the Cox proportional hazards model. This model contains a non-negative component, baseline hazard, which is not estimated or is poorly treated. In the present study, we aimed to determine survival prognostic factors in BC women using the new MPL method in the Cox model and compared it with the usual PL method to introduce the better method in identifying factors that really affect the survival of BC patients. MPL is a new method in which the regression coefficients and baseline hazard are estimated simultaneously using the new Newton-MI algorithm. This algorithm ensures that the estimated hazard function is non-negative. Based on our results, the standard deviations of most variables given by the new MPL method were less than those obtained by the usual PL method. These results are in line with the results of Ma et al. (11), Xu et al. (16) and Arena et al. (26). The penalty term in the MPL method helps estimate the baseline hazard more accurately and correct the bias of estimates (11); therefore, the final results and the hazard ratios will be more accurate.

In our study, in addition to lymphovascular invasion and recurrence variables, the estrogen receptor variable also played a significant role. The estrogen receptor is a predictor of survival in women with BC (27, 28) as confirmed by the new MPL method in our study. This is while this variable was not significant in the analysis using the usual PL method. Arena et al. used the MPL method to obtain more reliable variance parameter estimates (26).

Based on our findings, the new MPL method showed an association between lymphovascular invasion and survival as the risk of death was 2.46 times higher in patients with positive lymphovascular invasion than in patients with negative lymphovascular invasion (HR = 2.46). This outcome was in line with the report of Akbari et al. (27). One of the indices used to evaluate the quality of care is the estimation of the five-year survival rate. The five-year survival rate of BC patients in developed countries, regardless of the stage of cancer, is 73% (29). In this study, the MPL method calculated the overall five-year survival rate of 70% with a lower standard error than the PL method; this result is consistent with the results of Movahedi et al. (30) with a 71% survival rate.

In this study, the new MPL method calculated slightly smaller standard errors for most variables than the usual PL method and the estrogen receptor variable was a significant predictor of survival, as well. Although using the new MPL method does not guarantee lower standard errors for all estimates, since it estimates both baseline hazard function and coefficients simultaneously, the standard errors of most variables decrease. A simulation study by Xu et al. demonstrated that the MPL method works well and usually offers smaller standard errors and bias (15). Therefore, it can be applied to identify predictive factors more accurately than previous methods do.