Characterizing Trends in Cancer Patients' Survival Using the JPSurv Software

Abstract Background: Improvements in cancer survival are usually assessed by comparing survival in grouped years of diagnosis. To enhance analyses of survival trends, we present the joinpoint survival model webtool (JPSurv) that analyzes survival data by single year of diagnosis and estimates changes in survival trends and year-over-year trend measures. Methods: We apply JPSurv to relative survival data for individuals diagnosed with female breast cancer, melanoma cancer, non–Hodgkin lymphoma (NHL), and chronic myeloid leukemia (CML) between 1975 and 2015 in the Surveillance, Epidemiology, and End Results Program. We estimate the number and location of joinpoints and the trend measures and provide interpretation. Results: In general, relative survival has substantially improved at least since the mid-1990s for all cancer sites. The largest improvements in 5-year relative survival were observed for distant-stage melanoma after 2009, which increased by almost 3 survival percentage points for each subsequent year of diagnosis, followed by CML in 1995–2010, and NHL in 1995–2003. The modeling also showed that for patients diagnosed with CML after 1995 (compared with before), there was a greater decrease in the probability of dying of the disease in the 4th and 5th years after diagnosis compared with the initial years since diagnosis. Conclusions: The greatest increases in trends for distant melanoma, NHL, and CML coincided with the introduction of novel treatments, demonstrating the value of JPSurv for estimating and interpreting cancer survival trends. Impact: The JPSurv webtool provides a suite of estimates for analyzing trends in cancer survival that complement traditional descriptive survival analyses.


Introduction
During the past two decades, substantial progress has been made in the treatment, early detection and prevention of many cancers. Population-based cancer statistics, especially trends in cancer rates (incidence and mortality), are the measures used most frequently to report and monitor progress in cancer control (1,2). The Joinpoint regression model analyzes rates and proportions over time to identify timepoints (joinpoints) at which trends have changed, and to estimate the regression function with joinpoints identified (3). The Joinpoint regression model also provides a summary of the pace at which rates are changing, usually reported as the annual percent change (APC).
More recently, a Joinpoint survival model (4) was developed that allows for analysis of survival trends by single calendar year of diagnosis. The Joinpoint survival model (4) is an extension of the proportional hazards model for survival, where the effect of year of diagnosis is linear on the log of the probability of cancer death scale (5). Like the Joinpoint model for rates, the Joinpoint survival model estimates the location and number of calendar years where changes in survival trends occurred. However, the Joinpoint survival model has been underutilized in survival studies. Most studies investigating changes in cancer survival simply compare 5-year relative survival of patients with cancer diagnosed in earlier versus later grouped calendar years of diagnosis (6,7). This type of comparison does not provide information on when or how cancer survival is changing. Reasons for the limited use of the Joinpoint survival model may be related to challenges in model and trend measure interpretation as well as a lack of user-friendly software. Applying the Joinpoint survival model (4) requires knowledge of statistical software packages, for example, SAS or R. Importantly, the trend measure was previously defined as the percent change in the annual probability of cancer death (hazard) scale (4), which is a less used and less intuitive measure of a patient's prognosis.
To make the joinpoint survival analyses more readily available, we developed a JPSurv webtool that can be accessed at https:// analysistools.cancer.gov/jpsurv/. The JPSurv webtool is based on the JPSurv R package (https://cran.r-project.org/web/packages/ JPSurv/), which uses an iteratively reweighted least squares (IRLS) algorithm to estimate the Joinpoint survival model for relative survival, cause-specific and overall survival. To facilitate the communication of survival trends, we also developed a new measure summarizing trends on the survival scale.
We first present the JPSurv model and the trend measures on the hazard and survival scales. We then briefly describe the JPSurv webtool. We provide an application using relative survival data for patients diagnosed with female breast cancer, melanoma cancer, chronic myeloid leukemia (CML) and non-Hodgkin lymphoma (NHL). These cancer sites were chosen because they illustrate interpretation of survival changes and trends measures for a range of prognoses, and varying dissemination of screening and novel treatments.

Materials and Methods
The JPSurv software has been developed to analyze survival trends by single year of diagnosis (year-over-year). The Joinpoint survival model (4) is an extension of the proportional hazards model for survival, where the effect of calendar year of diagnosis is linear on the log hazard of cancer death scale. The JPSurv model allows for different linear trends between joinpoints.

The JPSurv model
We assume that the hazard rate of dying at a time t since diagnosis (follow-up) for people diagnosed at calendar year x, lðtjxÞ, follows a proportional hazards model where l 0 ðtÞ is the baseline hazard and hðxÞ is a joinpoint model, . . . ; d K the regression coefficients. The slope of the first segment before the joinpoint t 1 is b and the slope of the segment between joinpoints t k and t kþ1 is b þ P k j¼1 d j . When there is no joinpoint, then hðxÞ ¼ b x and the model becomes a Cox proportional hazards model with calendar year of diagnosis x as a covariate.

Projecting survival using the JPSurv model
The JPSurv model can also project survival beyond observed data. Thus, suppose data are calendar years 1 to L. Survival for patients diagnosed in calendar year x, x > L, can be estimated as SðtjxÞ ¼ S 0 ðtÞ expfhðxÞg , where S 0 ðtÞ ¼ expfÀ R t 0 l 0 ðuÞ dug and by extending the last joinpoint segment, that is, calculating The projection assumes that the trend in the last segment continues in future years. Note that for calendar year L usually only 1-year survival, Sð1jLÞ, is observed. Thus, the model also allows for estimation of 2-year, 3-year,. . ., 5-year survival for patients diagnosed at year L.

Model estimation and model selection
The JPSurv model uses a combination of grid search and IRLS methods (8) to estimate the number of joinpoints, the location of joinpoints, and the regression coefficients of the Joinpoint model. For each number of k joinpoints up to a maximum number of 5 joinpoints, we find the model (joinpoint location and regression parameters) that maximizes the likelihood function. Let M k denote the estimated kjoinpoint model and l k its maximum log-likelihood. The Bayesian Information Criterion (BIC) is defined as where n is the total number of follow-up years for all diagnosis years and p k is the number of parameters under the model. The Akaike Information Criterion (AIC) is also provided. Yu and colleagues (4) favored BIC and provided simulations and guidelines on the performance of BIC. In general BIC, had the best performance in selecting the correct model. However, when sample size was small it tended to detect a smaller number of joinpoints. The AIC is more sensitive to changes and estimates a larger number of joinpoints (4). SEs and confidence intervals can be calculated using the delta method as shown in Yu and colleagues (4). Details of the likelihood function and estimation methods for relative survival and cause-specific survival are provided in the Supplementary Materials and Methods.
Annual percent change in the probability of cancer death (APC_D) Let b Ã be the slope coefficient in a segment between consecutive joinpoints, so that hðxÞ ¼ b Ã x, then Yu and colleagues (4) suggested APC D ¼ 100 fexpðb Ã Þ À 1g% as the annual percent change in the hazard of cancer death for continuous follow-up time. For discrete data, cancer death hazard translate into l j ðxÞ the conditional probability of dying of cancer in interval j, being alive at the beginning of the interval, for people diagnosed at year x. The conditional probability of cancer death is also called interval or annual probability of cancer death.
Thus, an APC_D ¼ À2% represents an annual percent change of 2% decline in the annual probability of cancer death by single year of diagnosis. It represents a 0.98 annual risk of cancer death for patients diagnosed in a given year compared with those diagnosed in the prior year.

Annual absolute change in survival (AAC_S)
Prognosis is often reported as cumulative survival, that is, as the probability of being alive at 5 or 10 years after diagnosis. We developed a new measure to summarize survival trends on the cumulative survival scale. For calendar year x 2 ½t m ; t mþ1 between joinpoints with slope coefficient b Ã , the cumulative survival probability after t years from diagnosis (follow-up), is the product of the interval survival probabilities and can be written as expfÀl 0 ðsÞ expðb Ã xÞg: The absolute change in t-year survival at calendar year x is calculated as Because it varies by calendar year, we average over all calendar years in the interval. Thus, the (average) annual absolute change over a period of calendar years x i 2 ½t m ; t mþ1 Þ is estimated as ACSðt; x i Þ ðt mþ1 À t m Þ : The AAC_S(t) represents the average difference of cancer survival at t years from diagnosis for people diagnosed in a calendar year compared with people diagnosed in the prior year. Since it reflects the difference of 2 percentages, it is measured in percentage points. For example, moving up from 40% to 44% is a 4 percentage point increase in cumulative survival. Table 1 provides comparison and interpretation of APC_D and AAC_S.
The annual absolute change in survival, AAC SðtÞ, depends on follow-up survival time t. It can be calculated between any 2 calendar years to estimate the average year-to-year of diagnosis trend in the period, e.g., the most 10 recent calendar years. However, it is usually calculated in a segment between joinpoints. The SE for AAC SðtÞ is obtained using the delta method and is shown in the Supplementary Materials and Methods.

The JPSurv webtool
JPSurv is web-based tool available at (https://analysistools.cancer. gov/jpsurv/). Input data consist of grouped relative, cause-specific or overall survival data by annual time since diagnosis (follow-up) and by calendar year of diagnosis. Although more granular survival can be of interest for rapidly fatal cancer types, the tool requires annual survival rates and annual calendar year at diagnosis to calculated year-overyear changes in survival. The input data can be imported as a delimited text file or as a SEER Ã Stat text file and dictionary file calculated using the SEER Ã Stat software (https://seer.cancer.gov/seerstat/). The following variables are required: survival time interval in annual years since diagnosis, calendar year of diagnosis, number at risk at beginning of interval, number of cases lost to follow-up in the interval, number of cancer deaths (for cause-specific survival), number of deaths (for overall survival), and number of deaths and interval expected survival (for relative survival). Other stratifying variables such as cancer site, sex, stage, etc. which define survival cohorts can also be included. The user has full flexibility to define multiple cohorts by stratifying the survival calculations in SEER Ã Stat. For example, data may include survival for men and women diagnosed with colorectal and lung cancer. If the user does not select specific values for cancer site or sex, JPSurv will fit models to each combination of sex and cancer site in a single run using the same settings for all cohorts. If the user would like to change the settings for males lung cancer survival, than the user needs to select males and lung cancer, change the settings and run JPSurv for that cohort.
Settings include the maximum number of joinpoints tested (default is 0 joinpoints and maximum is 5), and the diagnosis years to be used in the model estimation (default is all years in the data). For two or more joinpoints or more than one survival cohort computational time is long and the user is required to submit an e-mail address to retrieve results. Advanced options provide control of the minimum number of years required between joinpoints, first and last intervals and the number of years for projections. The JPSurv help (https://analysistools.cancer.gov/jpsurv/html/help. html) and tutorial provide more detailed information on the input data and parameters to be specified by the user.
The default final model is the model with minimum BIC; however, users can use AIC to select their final model and all outputs can be displayed for any of the tested models. The outputs include three types of plots with predicted and observed values as below.
(i) Observed and modeled X-year relative survival (or cause-specific survival; y-axis) by calendar year of diagnosis (x-axis). The user can select multiple times since diagnosis X to display, for example, 1-year and 5-year. (ii) X to Xþ1 annual probabilities of cancer death (y-axis) by calendar of diagnosis (x-axis). The user can select multiple times since diagnosis X to display. (iii) Cumulative relative survival or cause-specific survival (y-axis) by time since diagnosis (x-axis) for people diagnosed in year Y. The user can select multiple calendar year of diagnosis and display for example, three survival curves by time since diagnosis for patients diagnosed in 1990, 1995, and 2000. The graphs allow for some customization. Trend measures are also calculated and can be displayed in the plots.

Application to the SEER data
We used data from the Surveillance, Epidemiology and End Results (SEER) program November 2019 data submission. We calculated relative survival for patients with cancer diagnosed between 1975-2015 and followed through 2017 in the SEER-9 areas using the SEER Ã Stat software (9). We selected patients diagnosed with female breast and melanoma cancers by historic stage (local, regional, and distant) at diagnosis (10) and NHL and CML. Because of changes in staging codes, historic stage is only available through cases diagnosed in 2015. We excluded cases diagnosed by autopsy or death certificate and those lost to follow-up in the month of diagnosis (alive and no survival time). Although patients diagnosed many years ago have more than 10 years of follow-up, we only Table 1. Comparison between measures to summarize trend: annual absolute change in survival (AAC_S) and the annual percent change in the conditional probability of cancer death (APC_D). Does it vary by time since diagnosis? Yes.

Motivation/Summary
More clear prognosis interpretation versus more awkward mathematical derivation. included relative survival for a maximum of 10 years of follow-up. The choice of follow-up time to include in the modeling depends on the interest in understanding longer survival. If only trends in 5-year survival will be reported, the user should only include up to 5-years of follow-up. We used the JPSurv software with the default parameters and tested a maximum number of 4 joinpoint for all the cancer sites and stage combinations with the exception NHL for which we used 5 joinpoints. The cancer sites in this study were chosen because of the introduction of new treatments that are known to have improved survival. We report models with the minimum BIC.

Clear mathematical interpretation versus challenging prognosis interpretation
Only for CML and to check the proportionality assumption we performed analyses limiting follow-up time to up to 2 and 5 years.

Overview of measures and interpretation
Relative survival has substantially improved since at least the mid-1990s for all combinations of cancer sites and stages (Fig. 1). The numbers near the lines represent the AAC_S between joinpoints for the 1-year, 5-year, and 10-year relative survival. Visually the largest improvements occurred for CML and NHL after the mid-1990s and distant melanoma after 2010. Supplementary Table S1 displays the BIC and AIC measures for all the fitted models including the chosen final model with the minimum BIC.
The trend measures and their respective 95% confidence intervals for the final models are displayed in Table 2 The trend measures in the survival (AAC_S) and in the annual probability of death (APC_D) scales have opposite signs ( Table 2). When AAC_S is positive and survival is increasing, then APC_D is   negative and the annual probability of dying of cancer is decreasing. In general, AAC_S and APC_D, provided consistent indication if the trends were statistically significant or not significant. No significant changes in survival trends or in the annual probabilities of cancer deaths were observed for local and regional breast cancer, local melanoma cancers in the late 1970s, for CML after 2010, or NHL in 1983-1990. The 1-year AAC_S was also non-significant for local-stage melanoma between 1981 and 2015. However, the AAC_S and APC_D measures differ in ranking the trends. The largest improvement in 5-year relative survival was observed for distant stage melanoma cancer after 2009, while the largest decline in the annual probability of cancer death was observed for local-stage breast cancer between 1982 and 1988. Because APC_D is a percent change (relative measure), it depends on the level of the probability of death. It reports the greatest change when the probability of death is small, and the prognosis is good. Figure 2 shows annual probabilities of cancer death and the respective APC_D measures on a variable y-axis scale. For local-and regional-stage breast and mela-noma cancer, the annual probabilities of dying of cancer are very small and it is challenging to plot using the 0 to 100 scale.

Model fit and the proportionality assumption: CML example
The joinpoint model did not fit well the CML survival data ( Fig. 1  and Fig. 2). The lack of fit, is more clearly understood from the annual probability of cancer death figures (Fig. 2). The probability of cancer death in the first year after diagnosis is underestimated and overestimated before and after 1995, respectively. The reverse occurs for patients who survive 4 years: the probability of cancer death in the 5th year after diagnosis is overestimated and underestimated before and after 1995, respectively. To check on the proportionality assumption, we performed analyses limiting follow-up time to up to 2 and 5 years. Using follow-up data up to 2 and 5 years, the BIC criteria identified respectively a model with 1 joinpoint in 1997 and 2 joinpoints in 1999 and 2010, (Table 3; Fig. 3). Figure 3 shows that the fit of the observed data to the 1-year annual probability of death between 2000-2010 is improved when using only 2 years of follow-up. The lack of fit using longer follow-up is because the baseline hazard of cancer death did not decrease proportionally before and after 1997-1999, which roughly coincides with the introduction of imatinib in 2000 to treat CML. In other words, the lack of proportionality shows that although the probability of dying of CML decreased for all intervals after diagnosis, it decreased more for patients surviving 3 or more years from diagnosis in the period after the introduction of imatinib compared with the period before. The AIC criteria identified 2 joinpoints in 1999 and 2010 using only 2 years of follow-up, which provides a joinpoint closer to the introduction of imatinib and better projections. Using 5 years the AIC and BIC select the same model ( Table 3).

JPSurv webtool and trend measures interpretation
We developed a user-friendly webtool that makes analyses of survival trends more broadly available without the need for any programming skills. We introduced a new measure on the cumulative survival scale, Table 3. Model fit for CML using up to 2 years and 5 years of follow-up data to fit the JPSurv model.

Implication of the SEER survival trend results
In applying JPSurv to the SEER data, we have considered cancer sites that were impacted by the introduction of novel treatment and screening (breast and melanoma). Survival increased for all cancer sites and all stage combinations. The largest improvement in survival occurred for distant-stage melanoma cancer after 2009, with an increase of 2.8 percentage points in 1-year relative survival for each subsequent year after diagnosis. This improvement coincides with the rapid dissemination of immune therapy for the treatment of advanced melanoma after 2011 (11)(12)(13)(14).
The second largest improvement in relative survival was observed for CML between 1995 and 2010. However, the Joinpoint survival model did not fit well the data. The observed survival data violated the baseline hazard proportionality assumption and did not decrease proportionally for all follow-up times across all years of diagnosis. Thus, the trend measure may not be accurate. The first tyrosine kinase inhibitor (TKI) to treat patients with CML was approved by the FDA in 2000 (15,16). Restricting the data to 2 years of follow-up, improved the fit of 1-year and 2-year survival and identified a joinpoint in 1997 using the BIC criteria and 2 joinpoints in 1999 (closer to 2000) and 2010 using the AIC criteria. The lack of fit indicated that after 2000 there were larger improvements in the chances of surviving CML for patients who had already survived 3 years or more compared with the initial years after diagnosis. This may indicate that patients with CML treated with TKIs who survive 2 or 3 years, subsequently experience a lower hazard of death and were in remission for a long time compared with patients diagnosed in the pre-TKI era.
Relative survival for NHL also increased substantially, 2.0 survival percentage points between 1995 and 2003, coinciding with the dissemination of CHOP (cyclophosphamide, doxorubicin, vincristine, and prednisone) therapy in 1993. This chemotherapy regimen, not only induced remissions in most patients with aggressive NHL, but it also "cured" a significant number of patients with some subtypes of NHL (17).

Challenges and opportunities
One important limitation of the JPSurv model is the proportionality assumption of the baseline hazard by time since diagnosis for different calendar years of diagnosis. For most applications and in all of the examples except CML, the proportionality assumption was reasonable (18). It assumes that a new treatment decreases the baseline hazard rate proportionally compared with a year in which conventional treatment was given. For CML this was not true. The observed relative risks of dying of cancer in the first year since diagnosis compared with the fifth year was on average 2 prior to the introduction of imatinib and became 4 after 2000, indicating a larger improvement for patients surviving 3 or more years. Relaxing the proportionality assumption adds complexity to the estimation and interpretation. The figure displaying the probability of death for different time since diagnosis can provide a visual check for the proportionality assumption. We are developing extensions of the JPSurv model that will relax the proportionality assumption. In the case of nonproportional data, the JPSurv model can still be applied to survival data with limited follow-up time and also provide insights into how treatment has differentially improved survival by time since diagnosis, as shown in the CML example.
The examples show that JPSurv can be helpful for understanding the impact of cancer treatment advances at the population level, especially for cancer sites unaffected by early detection or screening, such as: NHL, CML, and distant-stage melanoma. However, screening and early detection can introduce biases that artificially increase survival, which complicate interpretation. Screening can (i) advance the time of diagnosis (lead time bias); (ii) include a higher proportion of slower growing cancers, which are most likely to be picked up by screening (length bias); and (iii) detect slow growing cancers that would never cause symptoms or death (overdiagnosis). In the presence of screening, previous work (19,20) recommends analyzing changes in survival in conjunction with incidence and mortality trends. Thus, it is possible that the improvements in localized breast and melanoma cancers were attributable to screening biases rather than to real improvements in treatment.
The JPSurv model can also be used to predict survival beyond the observed data, both in future calendar years as well as time since diagnosis. Cancer patients' survival projections are useful to for estimating survival among more recently diagnosed patients, in simulation studies, and in prevalence projections to study sensitivity to different future survival projection assumptions. Previous work (21) has shown that extrapolating linear survival trends provides in most cases accurate survival predictions for up to 4 future years. However, caution should be used when extrapolating survival for longer number of years and when the last joinpoint is close to the last year of data. The BIC criteria, used as JPSurv default, is more parsimonious than the AIC and log-likelihood measures. However, users have flexibility to use less parsimonious models.
In summary, the JPSurv webtool provides a suite of estimates and graphs for analyzing cancer survival trends that complement traditional descriptive approaches. Our hope is that JPSurv will improve the reporting of cancer survival trends and interpretation. However, caution should still be used when interpreting survival trends for those cancer sites for which screening or early detection have been widely disseminated in the population, because survival increases may reflect biases rather than real improvements.