A Statistical Model with Non-Linear Effects and Non-Proportional Hazards for Breast Cancer Survival Analysis

Muditha Perera^*, Chris Tsokos

Department of Mathematics and Statistics, University of South Florida, Tampa, FL, USA.

DOI: 10.4236/abcr.2018.71005                

Abstract

The Cox proportional hazard model is being used extensively in oncology in studying the relationship between survival times and prognostic factors. The main question that needs to be addressed with respect to the applicability of the Cox PH model is whether the proportional hazard assumption is met. Failure to justify the subject assumption will lead to misleading results. In addition, identifying the correct functional form of the continuous covariates is an important aspect in the development of a Cox proportional hazard model. The purpose of this study is to develop an extended Cox regression model for breast cancer survival data which takes non-proportional hazards and non-linear effects that exist in prognostic factors into consideration. Non-proportional hazards and non-linear effects are detected using methods based on residuals. An extended Cox model with non-linear effects and time-varying effects is proposed to adjust the Cox proportional hazard model. Age and tumor size were found to have nonlinear effects. Progesterone receptor assay status and age violated the proportional hazard assumption in the Cox model. Quadratic effect of age and progesterone receptor assay status had hazard ratio that changes with time. We have introduced a statistical model to overcome the presence of the proportional hazard assumption violation for the Cox proportional hazard model for breast cancer data. The proposed extended model considers the time varying nature of the hazard ratio and non-linear effects of the covariates. Our improved Cox model gives a better insight on the hazard rates associated with the breast cancer risk factors.

Keywords

Breast Cancer, Cox Model, Non-Linear Effects, Non-Proportional Hazards

1. Introduction

According to the statistics currently cited by the American Cancer Society, about 252,710 women will receive a new diagnosis of breast cancer and about 40,610 women will die from breast cancer. Breast cancer is the second leading cause of cancer death in women. Breast cancer has been studied worldwide to improve survival by focusing on finding causes, reducing risks, developing new diagnostic tests and creating new treatment protocols [1] . Prognostic models play a major role in these studies. These models can be used to obtain estimates of risks of different adverse events such as death and recurrence based on the clinical, lifestyle, and socio-economic factors associated with the disease. This information acquired is important in designing treatment protocols, increasing disease awareness and preventing possible causes of cancer. A prognostic model will be a valuable tool only if it is developed carefully, evaluating the underlying model assumptions and inadequacies and determining if the most relevant model to address the study objectives is selected. The assessment of the adequacy of statistical models is only possible through the combination of several statistical analyses and proper investigation regarding the purposes for which the statistical model was initially conceptualized and developed for.

Cox proportional hazard (CPH) model [2] is a popular method that is being used in studying the relationship between survival times and explanatory variables. Even though there are many model adequacy methods that have been developed for the CPH model, usage of these methods does not seem to be very popular in applications of this model in real life data. The goal of this study is to develop an extended Cox regression model to breast cancer data by assessing and adjusting non-linear effects and non-proportional hazards that exist in the data. Different methods that can be used to assess non-proportionalities of the covariates are discussed. Also, the ways that can be used to assess the linearity of the continuous covariates and the approaches to identify the correct form of non-linear effects are presented.

2. Data

Female breast cancer patients of age 20 years and above who were diagnosed with invasive ductal carcinoma during the years 1990 to 2000 were extracted from the SEER breast cancer database for the present study [3] . This is the most common type of breast cancer and accounts for about 70% breast cancer incidence. The selected study data consists of a random sample of 1000 patients. Potential prognostic factors, including race and age of the patient, tumor size, lymph node status, extension of the tumor, tumor stage and outcome of progesterone receptor assay (PRA) were selected according to the current knowledge about the risks of cancer. Race of the patient is categorized to white, black and other. Age was measured at the diagnosis. Tumor size is the largest dimension or diameter of the primary tumor and it is measured at the diagnosis in millimeters. The variable lymph node status represents whether regional lymph nodes examined pathologically contain metastases. Lymph node-negative means the lymph nodes do not contain cancer and lymph node-positive means the lymph nodes contain cancer. Cancer cells in regional lymph nodes may mean that cancer is more likely to spread to other parts of the body. Stage variable in our study represents AJCC stage 3rd edition (1988-2003). It is known that some breast cancer cells need hormones to grow. These cancer cells have hormone receptors inside which are special proteins that when hormones attach to those, the cancer cells grow. A pathologist examines the cancer cells and determines whether they have many hormone receptors (hormone receptor-positive) or few or no hormone receptors (hormone receptor-negative). These hormones are estrogen and progesterone. Breast cancers that are estrogen-positive also tend to be progesterone positive, vice versa [4] . As our data showed a strong relationship between these two test results, we studied only progesterone receptor status (PRA). Survival time until cancer related death is the response variable of interest and death by other causes, lost to follow up or alive at the end of the recording period is considered as censored.

3. Methods

3.1. Cox Proportional Hazard Model

The Cox Proportional Hazard (CPH) model [2] is the most commonly used method of statistical modeling of survival data. It models the hazard of a subject at a given time with a given set of covariate values. Let t_i be the failure time for subject i, where $i=1,2,\cdots ,n$ . Then according to the CPH model, the hazard function for subject i at time t_i (>0) conditional on the set of covariates $Z_i=\left(Z_{1i},\cdots ,Z_{pi}\right)$ is given by

$h_i\left(t_i|z_i\right)=h_0\left(t\right)\exp \left({\beta }_1{z}_{1i}+\cdots +{\beta }_p{z}_{pi}\right)$ , (1)

where $h_0\left(t\right)$ is the baseline hazard function and denotes the hazard function when all covariate values take zero (reference values) and $\beta =\left({\beta }_1,\cdots ,{\beta }_p\right)$ are the corresponding regression coefficients for Z, the given covariates.

Let the model given in Equation (1) consist of one explanatory variable Z which takes values 1 (say, treatment) and 0 (say, control). Then the hazard ratio for a subject with covariate value 1 versus a subject with covariate value 0 at time t is given by

$HR\left(t\right)=\frac{h\left(t|Z=1\right)}{h\left(t|Z=0\right)}=\frac{h_0\left(t\right)\exp \left(\beta \right)}{h_0\left(t\right)}=\exp \left(\beta \right)$ . (2)

This implies that the ratio of the two hazards is a constant which does not depend on time, t. That is, the hazards of the two groups remain proportional over time. This is the key underlying assumption of the CPH model and is called the proportional hazards assumption.

3.2. Assessing the Model Adequacy

Once we fit a model to the data we need to verify whether the model adequately fits to the data before proceeding to the model interpretations. Failure to do so will result in incorrect decisions and conclusions about the data. Similar to standard regression models, the linearity assumption of continuous covariates and the existence of unusual/influential values need to be assessed in CPH model. Most importantly, the main underlying assumption of the CPH model needs to be checked for any violations. In some cases, the data will not satisfy the PH assumption and the use of the CPH model to describe the data without correcting the assumption violations will lead to misinterpreted results. Different types of residuals can be computed from the CPH model to cater to specific aspects of model adequacy, namely, Martingale, Schoenfeld and score residuals. Such residual based methods of assessing the overall model adequacy, correct functional form of the continuous covariates, unusual/influential observations and proportional hazards assumption are presented in the following sections. More details about the tests can be found in [5] and [6] .

3.2.1. Overall Goodness-of-Fit

Overall goodness of fit of the model can be assessed using Cox-Snell residual plot [7] . The idea is to plot Cox -Snell residuals versus the cumulative hazard function of the Cox-Snell residuals. The points on the plot should fall on a straight line with unit slope if the data fits the model well. However, the final decision of the model shouldn’t be taken solely on this plot. In practice, it has been found that Cox-Snell plot is not sensitive to small model inadequacies and not reliable in small sample sizes. Therefore, along with this overall goodness-of-fit test we should proceed to check separately for the situations where model inadequacies can occur in a CPH model.

3.2.2. Correct Functional Form

Identifying the correct functional form of the continuous covariates is a crucial step in model development. However, it is not practiced much in health data analysis. Cumulative Martingale residual plots with the continuous covariates are useful in assessing the linearity of the variables. The smoothed curve to the plot indicates whether the effect of the variable is linear or non-linear. In addition, this smoothed curve gives a hint on the functional form of the relationship of the covariate to the hazard. Another method that can be used to evaluate the linearity is to use the observed and simulated cumulative Martingale residuals [8] . Under the linearity assumption cumulative Martingale residuals can be approximated by zero mean Gaussian process. Hence, a comparison of observed and simulated cumulative Martingale residuals would reveal any departures from the linearity assumption of continuous covariate. The method is to use one thousand simulations of the cumulative Martingale residual paths and compute the proportion of times that the maximum absolute values of the simulated paths exceeds the maximum absolute value of the observed cumulative Martingale residual path. This value serves as the p-value for a supremum type of formal test linearity assumption. If the simulated paths exceed the observed path relatively few times, then it is an indication of the violation of the assumption.

3.2.3. Unusual and Influential Data Values

Detection of unusual data values and influential data values on the parameter estimates of the CPH model can be done using statistics similar to leverage and dfbeta in standard linear regression models. Score residuals have similar properties as leverage values. For continuous predictors, the further the value is from the mean, the larger the absolute value of the score residual is. Graphs of the score residuals and covariates aid in identifying any subjects with unusual data values. A statistic that is similar to dfbeta that approximately measures the difference between a particular coefficient value and the new coefficient if a value is removed from the sample can be computed for CPH model using score residuals and covariance matrix of the estimated coefficients [9] . This value is sometimes called scaled score residual and plots of these residuals and continuous covariates are useful to examine any subjects that influence the parameter estimates.

3.2.4. Proportional Hazard Assumption

Scaled Schoenfeld residuals [10] can be used to identify violations of proportional hazards assumption. The method is to include a coefficient that varies with time to the model instead of constant coefficient (β) in Equation (1). The time varying coefficient takes the form

${\beta }_j\left(t\right)={\beta }_j+{\gamma }_j{g}_j\left(t\right)$ (3)

where $g_j\left(t\right)$ is a function of time that the user has to specify. Approximated scaled Schoenfeld residuals have a mean at time t given by ${\gamma }_j{g}_j\left(t\right)$ . As a result, the plot of scaled Schoenfeld residuals versus time can be used to assess whether ${\gamma }_j$ zero is or not. That is, if slope is zero then $g_j\left(t\right)$ doesn’t depend on time and hence the hazard ratio is also constant with respect to time. In addition, a formal test to check whether ${\gamma }_j$ is zero has been proposed by [10] .

Another method that can be used is to use a transformation of Martingale residuals which is called score process [8] . Under the assumption of proportional hazards this process can be approximated by zero mean Gaussian process. Hence, a comparison of observed score process and simulated score processes under the PH assumption would reveal any departures from the assumption. The method is to use one thousand simulations of the score process and compute the proportion of times that the maximum absolute values of the simulated processes exceeds the maximum absolute value of the observed score process. This value serves as the p-value for a supremum type of formal test of PH assumption. If the simulated processes exceed the observed process relatively few times, then it is an indication of the violation of the assumption. In addition, graphs of these observed and simulated processes can be used to identify the departures from the proportional hazards.

3.3. Adjusting Non-Linear Effects of the Covariates

When continuous predictors are present, the common and convenient practice is to include them as categorical predictors or as linear predictors in the model being studied. Categorization of a continuous covariate might lead to subjective categorizations and loss of information. Also, if a continuous predictor is incorrectly included as a linear effect then it might lead to misleading conclusions from the model. When non-linear effects are detected, we should attempt to find the correct functional form of the effect or a function that closely follows the non-linear effect. In practice, most of the time non-linear effects do not take common parabolic nature. Therefore, more advanced transformations are needed to approximate the functional form of the covariates.

Fractional Polynomial Method

Fractional polynomials can be used to describe complex relationships between the outcome and continuous covariates. Authors of [11] have developed a model selection procedure to select the best fitting fractional polynomial for a given covariate. The method is to use one polynomial term model (FP1) and a two-term polynomial (FP2) to capture the pattern of the relationship between the covariate and the outcome. Then a deviance difference test is used to compare and choose the best fractional polynomial model. Default order of entering covariates to this procedure is based on the statistical significance with respect to p-value. Optionally, we can choose the order that variables enter. Also, we can specify certain continuous variables of interest in the model to be linear. Assume that we have one continuous covariate (X) which we aim to find the correct functional form and time variable T. Then, candidate fractional models can be written as

FP1 model: $h\left(x,t\right)=h_0\left(t\right)\exp \left\{{\beta }_1{x}^{p_1}\right\}$ (4)

FP2 model: $h\left(x,t\right)=h_0\left(t\right)\exp \left\{{\beta }_1{x}^{p_1}+{\beta }_2{x}^{p_2}\right\}$ (5)

where p₁ and p₂ are selected from the set (−2, −1, 0, 1, 2, 3); 0 corresponds to logarithm transformation. A summary of the process of fractional polynomial model building is given in Figure 1. More details about this method can be found in [11] .

3.4. Adjusting Non-Proportional Hazards Using Time Varying Effects Model

The common method that is used to account for non-proportionality in a covariate is stratification, that is, use of the proportional hazard violated variable as a grouping variable rather than a regressor in the model. Even though this method is simple and easy to understand it has some drawbacks. When stratified Cox model is fitted, it is not possible to estimate hazard ratios associated with the stratifying variable (non-PH variable). This will be a major limitation if the stratification variable is an important characteristic under the study. In addition, this method is more suitable for qualitative covariates as there will be loss of information if used for quantitative variable. Also, when the number of predictor variables that violates the proportional hazard assumption is large, stratified Cox model is not very useful.

Given the limitations of the stratified Cox model, we want to introduce an extended Cox model with time varying coefficients which can be used to address

Figure 1. Process of fractional polynomial model selection.

those limitations. The idea is to create a time dependent coefficient, $\beta \left(t\right)$ , for the covariate which violates the proportional hazard assumption. That is, $\beta \left(t\right)=\beta f\left(t\right)$ ; where $f\left(t\right)$ is a function of time to reflect the time varying nature of the hazard ratio under study. $f\left(t\right)$ could be based on the theoretical knowledge about the covariate or scaled Schoenfeld residuals with smoothed curves. The Cox model with time varying coefficients for i^th individual ( $i=1,2,\cdots ,n$ ) can be written in the form

$h_i\left(t\right)=h_o\left(t\right)\exp \left\{{\displaystyle {\sum }_{j=1}^p{\beta }_j\left(t\right)Z_{ij}}\right\}$ ; (6)

where $h_o\left(t\right)$ is the baseline hazard function, i.e. hazard function when all covariates (Z_j; _{$j=1,2,\cdots ,p$} ) takes the reference values at time = 0 (time at origin).

Recall that in the CPH model hazard ratio, $\frac{h_i\left(t\right)}{h_o\left(t\right)}$ can be obtained by $\exp (\beta j)$

which is constant over the time. In contrast, in the time varying coefficient Cox model, the hazard ratio is time dependent. That is, $\exp \left({\beta }_j\left(t\right)\right)$ is the relative hazard of two individuals at time t whose Z_j variable differs by one unit and the remaining variables take the same values for both of them.

4. Results

The mean and the standard deviation of follow up times of the patients are 10.5 years and 5.2 years respectively and median survival time is 11 years. Sixty five percent of the study sample were censored observations; that is, alive at the end of the follow up period, lost to follow up or death by other cause. Overall 5 years and 10 years survival probabilities are 80% and 70% respectively. Tumor size at diagnosis had a mean of 22 mm with a standard deviation of 18.4 mm and age at diagnosis had a mean of 58.2 years with a standard deviation of 13.5 years. Age at diagnosis was centered at the average for meaningful interpretations for the baseline survival probability. Initial evaluation of the covariates was done using the univariate CPH models and all the covariates were significant at 5% significance level. Next, a multivariate CPH model was developed using backward elimination process (removal significance level = 0.05). Only the variable extension was not significant in the model. This served as the selected standard CPH model where summary of the estimates are shown in Table 1.

The next step is to evaluate the adequacy of this standard CPH model in the aim of improving the model if there are any inadequacies present. First, overall model adequacy was assessed using Cox-Snell residuals. The corresponding Cox-Snell residual plot is shown in Figure 2 where the graph deviates from the reference line which goes through the origin. Since there is some evidence for overall model inadequacy, the next step was to explore what makes the model inadequate. Three main aspects of the model were assessed; namely linearity of the continuous covariates, unusual/influential data points and the proportional hazard assumption.

First, any unusual values and/or influential values were identified. Score residuals were computed for the standard CPH model and plotted against age and tumor size to identify whether there are any records that have values that deviate from the rest of the data to a great degree. Figure 3(a) and Figure 3(b) display the score residual plots for age at diagnosis and tumor size at diagnosis. It can be seen that there are two values far apart from the other values on the top right of the score residual plot for age. Also, there are four values that differ from the other values on the score residual plot for tumor size. Dfbeta and age and tumor size were plotted and shown in Figure 3(c) and Figure 3(d) to identify any strong influential values on the parameter estimates.

It appears that two data points in the plot for age and five data points on the plot for tumor size deviate from the rest of data points to a great extent. These identified values were further assessed to check what subjects correspond to these unusual behaviors and how they affect the parameter estimates. A model without these six identified extreme values was fitted and there was 53% reduction of the coefficient estimate for race-other term. Tumor size change was more

Table 1. Results of the standard cox proportional hazards model.

Figure 2. Cox-Snell residual plot for the standard cox proportional hazards model.

(a) (b) (c) (d)

Figure 3. Score residual plot for (a) Age; (b) Tumor size and dfbeta plots for (c) Age; (d) Tumor size.

than a 100% change (0.007 to 0.016) which is expected because five of the identified records had larger tumor sizes, greater than 130 mm. Breast tumor sizes are typically less than 50 mm and sometimes they can be more than 50 mm. However, observance of a tumor size that is greater than 130 mm clinically is possible but it is rare [12] . Also, some inconsistencies of the values for tumor size, lymph node status and stage can be found in these five data points. Hence, we decided to disregard the identified data points from the further analysis as they are unusual in the study data.

There are two continuous covariates that we aim to identify the correct functional form, namely age and tumor size at diagnosis. Martingale residuals for the null model without the predictors were computed and plotted with age and tumor size along with smoothed curve. Figure 4 and Figure 5 represent the corresponding smoothed residual plots for age and tumor size respectively. It is clear that age and tumor size have a non-linear relationship with estimated log hazard. Both covariates appear to have higher estimated log hazard as the covariate values increase. We further assessed the non-linear nature of these relationships in the aim of finding the best form of function that describes effects of age and tumor size on log hazard.

The method of fractional polynomials was used to capture the nature of the non-linear effects. Both age and tumor size revealed significant transformations which confirm the observation we obtained from Figure 4 and Figure 5. The

most appropriate transformation for age is $F{P}_{\text{age}}={\left(\frac{\text{agecentered}}{10}\right)}^2$ and for tumor size is $F{P}_{\text{size}}=\ln \left(\frac{\text{size}}{100}\right)$ .

The standard CPH models with linear terms and the fractional polynomial model with non-linear terms for age and tumor size were compared. Partial likelihood ratio test revealed test statistic of G = 3461.324 − 3380.917 = 80.407 with 2 degrees of freedom p-value of 2.68 × 10⁻¹⁸. Hence, the model with fractional polynomials for the non-linear effects is significantly different from the standard CPH models with linear terms.

Next, we used scaled Schoenfeld residual plots and simulated score residual plots to graphically assess the proportional hazards assumption of the CPH model adjusted for non-linear effects. Scaled Schoenfeld residual plots with smoothed curve are shown in Figure 6 for race-other, lymph node-positive, age² and PRA-positive. These plots indicate possible proportional hazard violations in the corresponding covariates. Age² seems to have an upward trend/non-linear. Almost all of the points in race-other plot lie around the horizontal line. There are some isolated points on top which might be the reason for the slight upward trend. Lymph node-unknown and stage III don’t seem to have a large trend or deviation from the horizontal line. The plot of PRA-positive shows a clear downward trend. It seems like an exponential decay and then leveling off as time increases. Plots for other covariates display any noticeable deviations from the horizontal line. That is, those covariates support proportional hazards

Figure 4. Smoothed martingale plot for age (smooth = 0.615).

assumption. We also performed the test of proportional hazards by [10] and the results are summarized in Table 2. Comparison of the observations from these plots with the test gives evidence for significant departures from proportional hazards for the variables, race-other, PRA-positive and age². Also, the test gives a significant result of proportional hazard assumption violation by stage III.

Observed and simulated score residuals plots for stage III, age², tumor size and PRA-positive are shown in Figure 7. Plots for other variables didn’t indicate evidence for violation of proportional hazard assumption. Observed paths for PRA-positive and tumor size clearly deviate from the cloud of simulated paths. Stage III also shows some slight deviations from the simulated paths. As we can observe only ten simulated paths in the graph compared to the thousand simulations done for each graph, it is difficult to make strong observations from these plots. Supremum tests of non-proportionality which consider all the simulated paths would give more accurate findings. The corresponding supremum test results are given in Table 3. These results confirms the non-proportional hazards observed under scaled Schoenfeld residuals and plots for PRA and age variable. In addition, this supremum test suggests that stage III and lymph node-unknown variables do not satisfy the proportional hazards assumption. Evidence for

Figure 5. Smoothed martingale plot for tumor size (smooth = 0.529).

proportional hazards assumption violation was clearly visible in the plots that we examined and in the tests that we performed for age² and PRA. Age² seemed to have a log hazard that has a linear upward trend and PRA-positive shows exponential type decay. Lymph node-positive, stage III and ln(tumor size) showed some evidence of non-proportionality.

In order to develop a model that accounts for these varying hazard ratios, we incorporated time varying coefficients to all these possible non-proportional covariates. We included an interaction term for coefficients of race-other, stage III, lymph node-positive and age² to vary with time linearly, $\beta \left(t\right)=\beta \times t$ . PRA was allowed to vary in time in two ways, continuously and discretely; we named them model A and model B respectively. The vector z includes all the covariates with time fixed effects and $\beta$ is the corresponding vector of model coefficients. Even though, ln(tumor size), race-other, lymph node-positive and stage III showed some weak evidence of non-proportional hazards, we included these terms in both models with time varying coefficients (Model A & Model B).

(a) (b) (c) (d)

Figure 6. Scaled schoenfeld residual plots.

However, time varying coefficients of these terms were not statistically significant in both models A and B.

Model A: Exponentially decaying effect for PRA where rate of decay (k) was estimated from the scaled Schoenfeld smoothed residual plots. ${\beta }_{\text{PRA}}\left(t\right)={\beta }_{\text{PRA}}{\text{e}}^{-kt}$ . The model takes the form

$\begin{array}{c}h\left(t,Z\left(t\right)\right)=h_0\left(t\right)\{\exp ({\beta }_{\text{age}}\left(t\right)\times {\text{age}}^2+{\beta }_{\text{PRA}\left(\text{positive}\right)}\left(t\right)\times \text{PRA}\left(\text{positive}\right)\\ +{\beta }_{\ln \left(\text{tumorsize}\right)}\left(t\right)\times \ln \left(\text{tumorsize}\right)+{\beta }_{\text{race}\left(\text{other}\right)}\left(t\right)\times \text{race}\left(\text{other}\right)\\ +{\beta }_{\text{lymphnode}\left(\text{positive}\right)}\left(t\right)\times \text{lymphnode}\left(\text{positive}\right)\\ +{\beta }_{\text{race}\left(\text{other}\right)}\left(t\right)\times \text{stageIII}+{\beta }^{\prime }z)\}\end{array}$

Model B: A model with piecewise hazard function for PRA was developed. That is, hazard ratio of PRA to vary discretely with time. Time scale was partitioned into 2 year intervals and five dummy variables were created to represent the piecewise effects of PRA as shown in Table 4.

Model B takes the form

Table 2. Test of proportional hazards by Grambsch & Therneau (1994).

Table 3. Test of Proportional Hazards by Lin et al. (1993).

Table 4. Dummy variables for PRA in Model B.

$\begin{array}{c}h\left(t,Z\left(t\right)\right)=h_0\left(t\right)\{\exp ({\beta }_{\text{age}}\left(t\right)\times {\text{age}}^2+{\beta }_{\text{PRA1}}\left(t\right)\times \text{PRA}\left(\text{positive}\right)\\ +{\beta }_{\text{PRA2}}\left(t\right)\times \text{PRA}\left(\text{positive}\right)+{\beta }_{\text{PRA3}}\left(t\right)\times \text{PRA}\left(\text{positive}\right)\\ +{\beta }_{\text{PRA4}}\left(t\right)\times \text{PRA}\left(\text{positive}\right)+{\beta }_{\ln \left(\text{tumorsize}\right)}\left(t\right)\times \ln \left(\text{tumorsize}\right)\\ +{\beta }_{\text{race}\left(\text{other}\right)}\left(t\right)\times \text{race}\left(\text{other}\right)+{\beta }_{\text{lymphnode}\left(\text{positive}\right)}\left(t\right)\\ \times \text{lymphnode}\left(\text{positive}\right)+{\beta }_{\text{race}\left(\text{other}\right)}\left(t\right)\times \text{stageIII}+{\beta }^{\prime }Z)\}\end{array}$

(a) (b) (c) (d)

Figure 7. Observed and simulated score residual paths.

Estimated hazard ratios for PRA from the Model B are given in Table 5. It can be seen that after 4 years hazard ratio for PRA-positive is approximately equal to one. That is, the estimated risk for PRA-positive and PRA-negative individuals is almost same after 4 years. By observing the p-values, it can be seen that hazard ratios for time intervals 0 - 2 years and 2 - 4 years are statistically significant. Hence, we decided to create three time intervals 0 - 2, 2 - 4 and >4 years instead of five intervals and refit the piecewise Cox model.

Modified Model B: New dummy variables for PRA are defined with two time interaction terms as below.

${\beta }_{\text{PRA1}}\left(t\right)=\{\begin{array}{l}1;0\le t<2\\ 0;t\ge 2\end{array}$

and

${\beta }_{\text{PRA2}}\left(t\right)=\{\begin{array}{l}1;2\le t<4\\ 0;0\le t<2\text{and}t\ge 4\end{array}$

Table 5. Estimated hazard ratios for PRA in model B.

Modified model B takes the form

$\begin{array}{c}h\left(t,Z\left(t\right)\right)=h_0\left(t\right)\{\exp ({\beta }_{\text{age}}\left(t\right)\times {\text{age}}_{\text{FP}}+{\beta }_{\text{PRA1}}\left(t\right)\times \text{PRA}\left(\text{positive}\right)\\ +{\beta }_{\text{PRA2}}\left(t\right)\times \text{PRA}\left(\text{positive}\right)+{\beta }_{\ln \left(\text{tumorsize}\right)}\left(t\right)\times \ln \left(\text{tumorsize}\right)\\ +{\beta }_{\text{race}\left(\text{other}\right)}\left(t\right)\times \text{race}\left(\text{other}\right)+{\beta }_{\text{lymphnode}\left(\text{positive}\right)}\left(t\right)\\ \times \text{lymphnode}\left(\text{positive}\right)+{\beta }_{\text{race}\left(\text{other}\right)}\left(t\right)\times \text{stageIII}+{\beta }^{\prime }Z)\}\end{array}$

We compared the AIC values of the two time varying coefficient models that we fitted where it was 3386.739 for model A and 3387.685 for modified model B. Therefore, according to the Akaike’s information criteria, both model fits are similar. Hence, both of these models were considered for our next step.

As the final step in model building we checked the significance of all possible two way interactions in both model A and modified model B. We added all the interaction terms to the non-linear and non-PH adjusted model and used backward elimination method to remove the nonsignificant interaction terms. Only $\text{PRA}\left(\text{postive}\right)\times \text{lymphnode}\left(\text{unknown}\right)$ term was left significant in Model A. To compare whether the interaction model makes a significant improvement than the main effects model, we performed likelihood ratio test. Partial likelihood ratio test statistic is $G=-2\left[-1681.37-\left(-1678.464\right)\right]=2.906$ with a p-value of 0.0886 from chi-square distribution with 1 degree of freedom.

Time varying piecewise model with interactions resulted significant interactions of lymph node-unknown with race-other, stage II, PRA-positive and tumor size and also interaction of age with race-other. To compare whether the interaction model makes a significant improvement than the main effects model, we performed likelihood ratio test. Partial likelihood ratio test statistic is $G=-2\left[-1681.37-\left(-1672.95\right)\right]=8.42$ with a p-value of 0.2970 from chi-square distribution with 5 degrees of freedom.

The improvements made by the interaction models are not significant at 5% significance level. Hence, considering law of parsimony, we decided to proceed with the models with non-linear terms and time varying coefficients but without covariate interactions.

A summary of parameter estimates for the two time varying models and the standard CPH model is given in Table A1. Estimated coefficients for race and lymph node are similar in all three models. For stage, the parameter estimates are higher in the standard CPH model than in the two extended Cox models. As seen from the log-likelihood and the AIC values, both extended models fit the data similarly resulting similar parameter estimates for all the linear and PH satisfied covariates.

For a continuous variable with a linear effect, hazard ratio is interpreted as the relative risk between two individuals with a unit difference of covariate values and it doesn’t depend on the actual values that the variables take. Hence, under the standard CPH model, the estimated hazard ratio for tumor size can be obtained as below. $\hat{\beta }z$ represents the linear predictor for the covariate values that are being held constant (at reference levels). To obtain a more interpretable hazard ratio, let’s consider a 10 unit difference in tumor size. This will result in a hazard ratio of

$\begin{array}{c}\frac{\widehat{h\left(t,\text{size}=20\right)}}{\widehat{h\left(t,\text{size}=10\right)}}=\frac{h_0\left(t\right)\times \exp \left(0.0069\times 20+\hat{\beta }z\right)}{h_0\left(t\right)\times \exp \left(0.0069\times 10+\hat{\beta }z\right)}\\ =\exp \left(0.6176\times \left(20-10\right)\right)\\ =1.07\end{array}$

which indicates a 7% increase in hazard. Under the standard CPH model, this hazard ratio stays the same for any 10 unit increase in tumor size.

In contrast, when there is a non-linear effect present, hazard ratio depends on the covariate values that we are interested in and not only on the difference between values. How the hazard ratio is computed for the non-linear effect of tumor size under the piecewise Cox model is given below.

The parameter estimate for the ln(tumor size) from the piecewise Cox model is 0.6176. Say we need to estimate the hazard ratio between two individuals with tumor size 10 mm and 20 mm given that other covariate values are the same for both of them. Then the hazard function at tumor size = 10 mm is,

$\widehat{h\left(t,\text{size}=10\right)}=h_0\left(t\right)\times \exp \left(0.6176\times \ln \left(10\right)+\hat{\beta }z\right)$

and for tumor size = 20 mm,

$\widehat{h\left(t,\text{size}=20\right)}=h_0\left(t\right)\times \exp \left(0.6176\times \ln \left(20\right)+\hat{\beta }z\right)$

Now, the ratio of the two hazards is

$\begin{array}{c}\frac{\widehat{h\left(t,\text{size}=20\right)}}{\widehat{h\left(t,\text{size}=10\right)}}=\frac{h_0\left(t\right)\times \exp \left(0.6176\times \ln \left(20\right)+\hat{\beta }z\right)}{h_0\left(t\right)\times \exp \left(0.6176\times \ln \left(10\right)+\hat{\beta }z\right)}\\ =\exp \left(0.6176\times \ln \left(\frac{20}{10}\right)\right)=1.53\end{array}$

That is, for a 10 unit increase in tumor size from 10 mm, the risk of cancer death increases by about 50%. Unlike the hazard ratio estimated by the standard CPH model, the non-linear effect of tumor size results different hazard ratios for different tumor sizes that are being compared. It can be shown that hazard ratio for a 10 unit increase in tumor size from 20 mm to 30 mm is 1.28 and from 30 mm to 40 mm, it is 1.19.

In the standard CPH model where we assumed PH satisfied for PRA, the hazard ratio for an individual with PRA-positive compared to an individual with PRA-negative is $\exp \left(0.4303\right)=1.53$ . However, we found evidence that PRA-positive violated the PH assumption. That is, it doesn’t have a hazard ratio constant over time. Also, according to the literature on PRA status of women with breast cancer, we can find clinical evidence for this observation [13] .

In model A, we fitted a continuous function of time for the effect of PRA. That means at each point of time it results a different hazard ratio for PRA-positive relative to PRA-negative.

According to model A, estimated hazard ratio at time t is given by

${\text{HR}}_{\text{PRA}\left(\text{positive}\right)}=\exp \left(1.3726\times \exp \left(-0.23\ast t\right)\right)$

In modified model B, we modeled the effect of PRA in a piecewise time varying manner. For each 2 year interval from the start, we let the model estimate different coefficient for PRA-positive. After 4 years the hazards for PRA-positive and PRA-negative was not significantly different where it was near 1. Using the parameter estimates from Table A1, the hazard ratio for PRA-positive at time t is given by

${\text{HR}}_{\text{PRA}\left(\text{positive}\right)}=\{\begin{array}{l}\exp \left(0.0543+1.3152\right);0\le t<2\\ \exp \left(0.0543+0.7182\right);2\le t<4\\ \exp \left(0.0543\right);4\le t\end{array}$

Table 6 presents the estimated time varying hazard ratios for PRA-positive. At the start time, both models A and modified B seem to estimate the relative risk similarly with hazard ratios 3.95 and 3.93 respectively. Overall, model A estimates are higher than the modified model B estimates. Piecewise Cox model approaches to HR = 1 faster than the continuous time varying Cox model. Initially, an individual with PRA-positive has about four times of risk of cancer death than an individual with PRA-negative. At time equal to 2 years, the risk of cancer death for a PRA-positive individual is about twice of a PRA-negative person.

Because age has a non-linear effect and non-PH effect, special attention should be given when obtaining hazard ratios for age. We present the corresponding computations below for the piecewise Cox model. Say we need to estimate the hazard ratio between a 68.2 years old individual and 58.2 years old (mean age) individual given that other covariate values are the same for both of them. Then the hazard function at age = 58.2 years is

$\begin{array}{l}\widehat{h\left(t,\text{age}=58.2\right)}\\ =h_0\left(t\right)\exp \left\{0.0422\times {\left(\frac{58.2-58.2}{10}\right)}^2+0.0038{\left(\frac{58.2-58.2}{10}\right)}^2\times t+\hat{\beta }z\right\}\\ =h_0\left(t\right)\exp \left\{\hat{\beta }z\right\}\end{array}$

Table 6. Estimated time-varying hazard ratios for PRA-positive.

and at age = 68.2 years is

$\begin{array}{l}\widehat{h\left(t,\text{age}=68.2\right)}\\ =h_0\left(t\right)\exp \left\{0.0422\times {\left(\frac{68.2-58.2}{10}\right)}^2+0.0038{\left(\frac{68.2-58.2}{10}\right)}^2\times t+\hat{\beta }z\right\}\end{array}$

Now, the ratio of the two hazards is

$\begin{array}{l}\frac{\widehat{h\left(t,\text{age}=68.2\right)}}{\widehat{h\left(t,\text{age}=58.2\right)}}\\ =\frac{h_0\left(t\right)\exp \left\{0.0422\times {\left(\frac{68.2-58.2}{10}\right)}^2+0.0038{\left(\frac{68.2-58.2}{10}\right)}^2\times t+\hat{\beta }z\right\}}{h_0\left(t\right)\exp \left\{\hat{\beta }z\right\}}\\ =\exp \left\{0.0422\times {\left(\frac{68.2-58.2}{10}\right)}^2+0.0038{\left(\frac{68.2-58.2}{10}\right)}^2\times t\right\}\\ =\exp \left\{0.0422+0.0038\times t\right\}\end{array}$

Therefore, we get a time dependent expression for the relative hazard for two individuals for a 10 year increase in age from the mean age. Using this expression we computed relative hazards for different times and the results are given in Table 7. It can be seen that as time increases the hazard ratios are increasing at a slower rate. However, under the standard CPH model, relative risk for an increase of 10 years from mean age is $\exp \left(0.0381\times 10\right)=1.46$ irrespective of the time. When the model is adjusted for non-linear and non-PH effects we get different risk ratios than we would get from unadjusted Cox model.

To further visualize how the adjustments to the Cox model make the hazard estimations different, we graph hazard ratios with respect to the age increments and time as shown in Figure 8. It shows that at time = 0 hazard ratios for age are increased with age. However, approximately from time = 5 years, lower ages have higher risk than the mean age 58.2 years (baseline). When age is higher than the mean age, risks are increasing rapidly.

Figure 8 clearly shows how hazard ratios for age changes linearly with time and quadratically with age. Therefore, if we had used the standard CPH model to estimate the hazard ratios for age, it would not provide a flexible hazard ratio

Table 7. Estimated hazard ratios for age under the Cox model with piecewise time varying effects (Modified model B).

Figure 8. Hazard ratio plot for age adjusted for non-linearity and non-proportionality.

function as our extended Cox model which could explain the risk of cancer deaths more closely to the true pattern.

Recall that the summary of the estimated standard and extended Cox model parameters is given in Table A1. The following interpretations can be obtained by those estimates.

Under the piecewise Cox model (modified model B), risk of cancer death of a subject in race-black is about exp(0.5179) = 1.7 times higher than a subject in race-white. In contrast, hazard for a subject in race-other compared to a subject in race-white is 0.5. That is, race-white breast cancer patient is two times likely to have a death from cancer than a patient in race-other. Risk of cancer death for breast cancer patients gets larger as the stage of the disease gets higher as seen from Table A1 where hazard ratios for stage II, III and IV relative to stage I are 1.3, 1.8 and 6.2 respectively. All these estimates for the time-fixed variables are close to the estimates of the standard CPH model except for stage II and III. Also, these estimates are very close to the hazard ratio estimates from model A.

5. Discussion

Evaluating the model assumptions in CPH model should be an integral part of model building. Otherwise, it will result in incorrect conclusions about the data. Once model inadequacies are identified, they should be addressed before making any model interpretations. One of the important aspects to focus on when building a CPH model is the correct functional form of the continuous covariates. The simplest method is to categorize these covariates and include them in the model. This will lead to loss of information present in the covariate and the model would not provide accurate estimations about the survival probabilities. Including a continuous covariate as a linear effect is better than including it as a categorical variable. However, if the true effect is non-linear, then identifying the correct form of the effect will increase the predictive accuracy of the model. Fractional polynomial method [11] that we used in the present study is a good method to select the most appropriate functional form for the continuous covariates. Also, spline functions, which are piecewise polynomials connected across intervals of a given continuous covariate, are useful in finding the underlying non-linear effects of the covariates [14] .

Standard Cox model provide hazard ratios that are constant over the time. However, if there are covariates which have hazard ratios that change with time, then the standard CPH model will misinterpret the data. Hence, assessing the proportional hazard assumption in CPH model will not only help to decide the validity of the model, but also will provide a better understanding about the covariate effects. The major assumption of the CPH model, proportional hazards, is not easy to evaluate correctly. In some situations, the presence of other model inadequacies such as influential values and non-linear effects may cause proportional hazards tests to reveal significant non-proportionalities when actually they are proportional. Therefore, one first should assess and adjust the model for other inadequacies before performing proportional hazards tests.

Non-proportional hazards are more evident when longer follow up times are studied. If one is not interested in studying time effects, considering a shorter follow-up time might avoid the issue of dealing with non-proportional hazards. A simple and easy method to take non-proportional hazards into account is to use stratified Cox model. That is, stratify the CPH model by the non-proportional covariate. However, this method is mostly suitable when there is one non-PH variable and it is a categorical variable. Also, when the CPH model is stratified by a covariate, we are unable to estimate the corresponding effect on the survival probability. A more flexible method to handle non-proportional hazards is the use of extended Cox regression model with time varying effects/hazard ratios. This model can address more than one non-PH variable both continuous and categorical. Also, it provides estimates of time varying hazard ratios. This is the approach we used in this study to address non-proportional hazards. One difficulty of using this model is to find the correct function of time to include in the time varying coefficients. Schoenfeld residual plots used to detect proportional hazard violations might guide to find this function of time. Depending on the study objectives and clinical plausibility time varying effects can be modeled as piecewise functions or continuous functions.

In the present study, we identified and adjusted non-linearities and non-proportionalities that exist in the breast cancer data. The model development started by fitting the standard Cox model and then checked for the model inadequacies: non-linear effects, non-proportional hazards and unusual/influential values. Our analysis suggested few unusual and influential data points. We performed a sensitivity analysis of the parameter estimates with and without these points and found out that removal of these points changes the estimates of race-other and tumor size by more than 50%. Also, there were some inconsistent values taken by lymph node, extension and stage variables among the data points in this identified list. Therefore, these data points were not considered for further analysis of the current study. Our model building process revealed non-linear effects in both of the continuous covariates that we considered. The method of fractional polynomials proposed a logarithm effect for tumor size at diagnosis and quadratic effect for age at diagnosis. Our finding of a quadratic effect was consistent with findings of a similar study of breast cancer [15] . This effect suggests higher risk of cancer death for younger females and older female relative to middle aged females.

PRA and age were found to be violating the proportional hazards assumption under all the evaluation methods that were considered. Non-proportionality of the PRA was modelled through a continuous time dependent function guided by scaled Schoenfeld residual plot. Also, a piecewise time dependent function was used to model the time dependency of the PRA effect. In both of these models the effect of age was modelled through a time dependent quadratic effect. Both of these extended Cox models had very close log-likelihood values and AIC values. However, the estimated hazard ratios were fairly different for PRA under these two competing models. Up to 4 years both models gave similar relative risks of We found that under the piecewise effect Cox model that risk of PRA-positive relative to PRA-negative is not significant after around 4 years. This finding is consistent with the results of similar studies on breast cancer [13] where they discuss that the difference of the effects of PRA-positive and PRA-negative diminishes after around 5 years. According to our continuous time varying effect model, it seemed the differences in the risk decrease but at a lower rate and it approaches 1 approximately 13 years from the date of diagnosis. Considering this fact, the piecewise Cox model is more suitable for the data being studied. The effect that we found for age is interesting in that it contained both non-linear and non-proportional hazards. According to our extended Cox model, age had a linear effect on hazard ratio up to around 3 years and after that it shows a quadratic effect (Figure 8). For example, given that the time is more than 6 years from the diagnosis, an individual as young as 28.2 years old has at least the same risk of cancer death as a 68.2 years old individual. This result is consistent with a similar discussion where they suggested a higher risk of breast cancer for younger and older female than the middle aged females [15] .

6. Conclusion

We have identified that effects of age and tumor size at diagnosis on the hazard function are quadratic and logarithmic respectively. Also, we found that age and PRA-positive violate the assumption of proportional hazards. To address all these inadequacies of the standard Cox model, we have developed a more flexible extended Cox model with non-linear effects for age and tumor size and with non-PH effects of PRA and age described by a piecewise time dependent coefficient and linear time dependent coefficient respectively. This model gives improved and more accurate estimates of the risks of cancer specific death for women diagnosed with breast cancer.

Appendix

Table A1. A Comparison of initial and the extended Cox proportional hazards models on breast cancer data.

FP_size = ln(tumor size) and FP_age = age².

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	American Cancer Society (2016) Breast Cancer. https://www.cancer.org/cancer/breast-cancer/about/whats-new-in-breast-cancer-research.html
[2]	Cox, D.R. (1972) Regression Models and Life Tables. Journal of Royal Statistical Society, 34, 187-220.
[3]	National Cancer Institute (2012) Surveillance, Epidemiology, and End Results (SEER) Program. Research Data (1973-2009). National Cancer Institute, Surveillance Research Program, Surveillance Systems Branch. http://www.seer.cancer.gov
[4]	Joslyn, S., Gesme, D. and Lynch, C. (1996) Estrogen and Progesterone Receptors in Primary Breast Cancer. The Breast Journal, 2, 187-196. https://doi.org/10.1111/j.1524-4741.1996.tb00094.x
[5]	Collett, D. (2003) Modelling Survival Data in Medical Research. Chapman Hall, Boca Raton.
[6]	Hosmer, D., Lemeshow, S. and May, S. (2008) Applied Survival Analysis: Regression Modeling of Time to Event Data. 2nd Edition, John and Wiley Sons Inc., Hoboken. https://doi.org/10.1002/9780470258019
[7]	Cox, D. and Snell, J. (1986) Applied Statistics. Chapman and Hall, London.
[8]	Lin, D., Wei, L. and Ying, Z. (1993) Checking the Cox Model with Cumulative Sums of Martingale-Based Residuals. Biometrika, 80, 557-572. https://doi.org/10.1093/biomet/80.3.557
[9]	Cain, K. and Lange, N. (1984) Approximate Case Influence for the Proportional Hazards Regression Model with Censored Data. Biometrics, 40, 493-499. https://doi.org/10.2307/2531402
[10]	Grambsch, P. and Therneau, T. (1994) Proportional Hazards Tests and Diagnostics Based on Weighted Residuals. Biometrika, 81, 515-526. https://doi.org/10.1093/biomet/81.3.515
[11]	Sauerbrei, W. and Royston, P. (1999) Building Multivariable Prognostic and Diagnostic Models: Transformation of the Predictors by Fractional Polynomials. Journal of the Royal Statistical Society, Series A, 162, 71-94. https://doi.org/10.1111/1467-985X.00122
[12]	Breast Cancer Conditions. http://www.mayoclinic.org/diseases-conditions/breast-cancer/multimedia/tumor-size/img-20006260
[13]	Harris, J., Lippman, M. and Osborne, C. (2014) Diseases of the Breast. Lippincott Williams & Wilkins, Philadelphia.
[14]	Stone, C. and Koo, C. (1985) Additive Splines in Statistics. American Statistical Association, Washington DC.
[15]	Wingo, P., Ries, L., Parker, S. and Heath, C. (1998) Long-Term Cancer Patient Survival in the United States. Cancer Epidemiology, Biomarkers & Prevention, 7, 271-282.