Baade et al. BMC Cancer (2015) 15:27 DOI 10.1186/s12885-015-1024-4 RESEARCH ARTICLE Open Access Prognostic survival model for people diagnosed with invasive cutaneous melanoma Peter D Baade1,2,3*, Patrick Royston4, Philipa H Youl1,2,3, Martin A Weinstock5,6,7, Alan Geller8 and Joanne F Aitken1,2,3 Abstract Background: The ability of medical practitioners to communicate risk estimates effectively to patients diagnosed with melanoma relies on accurate information about prognostic factors and their impact on survival. This study reports the development of one of the few melanoma prognostic models, called the Melanoma Severity Index (MSI), based on population-based cancer registry data. Methods: Data from the Queensland Cancer Registry for people (20–89 years) diagnosed with a single invasive melanoma between 1995 and 2008 (n = 28,654; 1,700 melanoma deaths). Additional clinical information about metastasis, ulceration and positive lymph nodes was manually extracted from pathology forms. Flexible parametric survival models were combined with multivariable fractional polynomial for selecting variables and transformations of continuous variables. Multiple imputation was used for missing covariate values. Results: The MSI contained the variables thickness (transformed, explained 40.6% of variation in survival), body site (additional 1.9% in variation), metastasis (1.8%), positive nodes (0.7%), ulceration (1.3%), age (1.1%). Royston and Sauerbrei’s D statistic (measure of discrimination) was 1.50 (95% CI = 1.44, 1.56) and the corresponding RD2 (measure of explained variation) was 0.47 (0.45, 0.49), demonstrating strong explanatory performance. The Harrell-C statistic was 0.88 (0.88, 0.89). Lacking an external validation dataset, we applied internal-external cross validation to demonstrate the consistency of the prognostic information across geographically-defined subsets of the cohort. Conclusions: The MSI provides good ability to predict survival for melanoma patients. Beyond the immediate clinical use, the MSI may have important public health and research applications for evaluations of public health interventions aimed at reducing deaths from melanoma. Keywords: Melanoma, Survival, Prognostic model, Thickness, Population-based, Risk Background The incidence of cutaneous melanoma in Australia is the highest in the world due to the combination of high ultraviolet radiation, outdoor lifestyle and predominately Caucasian population [1]. As most melanomas in Australia are diagnosed when thin, [2,3] overall survival from melanoma is high. Internationally, overall five-year survival estimates in Western countries exceed 85% [3-6], although they are lower in Eastern Europe [7]. However * Correspondence: peterbaade@cancerqld.org.au 1 Cancer Council Queensland, 553 Gregory Terrace, Fortitude Valley, Queensland, Australia 2 School of Public Health and Social Work, Queensland University of Technology, Brisbane, Queensland, Australia Full list of author information is available at the end of the article there remains an important subset of melanomas that are diagnosed at an advanced stage, many with clinically apparent metastatic spread, for which the prognosis is poor [3,8,9]. Prognostic models are used to assist clinicians with their assessment of a patient’s future outcome and so enhance their informed decision making process with the patient [10]. Most prognostic models for melanoma are derived from selected clinical cohorts composed of patients referred to tertiary hospitals or specialist cancer centres [11,12] and/or focussed on specific subgroups of melanomas such as localised [13] or nodular melanoma [14]. A similar limitation in regards to specific cohorts of melanoma patients holds for population-based studies of survival outcomes [2,15,16]. © 2015 Baade et al.; licensee BioMed Central. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Baade et al. BMC Cancer (2015) 15:27 Page 2 of 13 We describe here the development of a prognostic model for deaths from melanoma among all patients diagnosed with a single invasive cutaneous melanoma in Queensland, Australia, using data from a populationbased cancer registry, and discuss the potential application in clinical practice and research. performed using Stata/SE version 12.1 for Windows (StataCorp, TX, USA). The Royston-Parmar models were fitted using the stpm2 package [20,21]. Multiple imputation by chained equations Methods Ethics approval to conduct this study was obtained from the University of Queensland Behavioural & Social Sciences Ethical Review Committee. Since these data are not publically available, approval to access the required data was obtained from the Queensland Department of Health. Data were obtained from the Queensland Cancer Registry (QCR) for all patients in Queensland with a histologically confirmed diagnosis of invasive melanoma (C44, M872-879) for the period 1995 to 2008. Notification of all cancers, apart from squamous and basal cell cancer, is required by law. Melanomas diagnosed on the basis of metastasis only were excluded. Variables extracted for each patient were sex, age at diagnosis, date of diagnosis, anatomic sub-site of the melanoma, tumour thickness, (Clarks) level of invasion and tumour morphology. Information on ulceration, presence and extent of metastasis at diagnosis and the number of positive lymph nodes (local, regional or distal)at diagnosis was extracted from pathology forms held by the Registry. Due to the difficulty of attributing death to a specific melanoma, we excluded all patients who were known to have had more than one histologically confirmed diagnosis of invasive melanoma since the establishment of the QCR in 1982. As with a previous report [2] only patients aged 15–89 years at diagnosis were initially considered as there is evidence that melanoma survival outcomes are different for younger age groups [17] and death certificates are less precise for older patients [18]. The QCR database was matched against the Queensland Register of Births, Deaths and Marriages and the National Death Index to identify deaths in Queensland and interstate respectively up to 31st December 2010. Cause of death was coded using all available information from death certificates, autopsy reports and pathology reports. Melanoma-specific survival was estimated from the mortality of people diagnosed with melanoma between 1995 and 2008 (inclusive), with follow-up for all cases to December 31, 2010. Patients who died before 31st December 2010 from conditions other than melanoma were censored at date of death. Those who were not recorded as dying by 31st December 2010 were censored at this date. Median follow-up time was calculated using the reverse Kaplan-Meier method [19]. All data analyses were Due to missing values in the included variables (Table 1), a complete case analysis would have excluded at least 30% of the initial cohort, potentially introducing a bias if the excluded cases were a non-random sample. We used multiple imputation [22] methods to deal with the missing data, using the mi impute chained and mi estimate commands for chained equations and subsequent regression model estimation. In the imputation modelling we included the Kaplan-Meier estimate of the survival curve, vital status, and the interaction between the survival and vital status. All variables included in the prognostic model were also included in the series of chained imputation equations. We used n = 30 imputations based on the percentage of incomplete cases [23]. Predictive mean matching was used for the imputation of thickness, logit models for lymph nodes (none vs one or more) and ulceration (yes, no), and multinomial logit models for Clark’s level and subsite. Derivation of the survival model Kaplan-Meier survival estimates, stratified by each of the covariates, were calculated to describe the melanoma cohort. We then developed a multivariable survival model to generate the prognostic index. Multivariable survival analyses are generally carried out using Cox regression. However several authors have highlighted the limitations of this method for prognostic models, [21,24] particularly relating to the appropriate modelling of the baseline hazards function. A parametric alternative to the Cox model, known as a flexible parametric survival model, is fitted on the log cumulative hazard scale [20,21]. These Royston-Parmar (RP) models use natural cubic splines to estimate the baseline cumulative hazard function. The selection of scales and number of degrees of freedom for the baseline spline function was made based on the Bayes information criterion (BIC) statistic. For our data, the probit scale with 3 degrees of freedom provided the best fit. These 3 degrees of freedom equate to 2 interior knots along with the 2 boundary knots. We then ran a multivariable fractional polynomial (MFP) procedure designed for multiple imputed data to help guide the backwards selection of covariates, and identify appropriate transformation(s), to retain in the fully adjusted survival model [25,26]. To improve the fit, we incorporated an additional smooth rank transformation [27] which accommodates sigmoid dose–response relationships. Visual inspection of the smoothed martingale residuals for the regression models was also used to Baade et al. BMC Cancer (2015) 15:27 Page 3 of 13 Table 1 Cause-specific survival: invasive Melanoma in Queensland (1995–2008, follow-up to 2010) Total Total Gender Male Female Age group (years)a 20-39 40-59 60-69 70-79 80-89 Thicknessa <=0.50 mm 0.51-1.00 1.01-1.50 1.51-2.00 2.01-4.00 4.01+ Unknown Body site Scalp Face/Lip/Eyelid Ear/Neck Chest/Axilla Abdomen + Hip Back/Buttocks Upper arm/Forearm/Hand/Finger/sb.hand Foot/sb. foot/heel/toec Thigh/leg/ankle Shoulder Unknown b b 1 year survival 99.1 5 year survival 95.0 10 year survival 92.6 28,654 (100%) 16,269 (57%) 12,385 (43%) 98.9 99.4 93.8 96.5 90.8 94.8 4,958 (17%) 10,814 (38%) 5,405 (19%) 5,006 (17%) 2,471 (9%) 99.8 99.5 99.1 98.7 97.5 97.4 96.5 94.6 92.0 88.6 96.1 94.4 91.7 87.9 85.2 12,641 (44%) 8,107 (28%) 2,487 (9%) 1,386 (5%) 2,234 (8%) 1,170 (4%) 629 (2%) 100 99.9 99.2 98.8 97.3 92.9 91.5 99.6 97.9 93.1 87.9 79.6 67.8 82.6 99.0 96.2 88.3 80.3 70.8 61.4 78.6 565 (2%) 2,221 (8%) 1,741 (6%) 1,435 (5%) 791 (3%) 8,046 (28%) 3,644 (13%) 384 (1%) 5,961 (21%) 2,452 (9%) 1,414 (5%) 96.2 99.0 98.9 99.2 98.7 99.5 99.5 97.4 99.5 99.5 96.4 78.8 94.5 93.4 96.0 93.7 95.2 96.6 87.7 96.5 96.0 92.0 69.3 90.4 90.6 94.0 90.7 92.8 95.4 79.6 94.7 93.7 90.2 Morphology Superficial spreading melanoma Nodular melanoma Malignant melanoma in junctional naevus Lentigo Maligna melanoma Acral lentiginous melanoma Other specified melanoma Malignant melanoma NOS Ulceration No ulceration Known ulceration No mention on pathologyb 17,110 (60%) 2,627 (9%) 8,917 (31%) 99.7 95.4 99.1 97.3 77.0 95.5 95.5 70.7 93.1 16,229 (57%) 2,415 (8%) 654 (2%) 1,688 (6%) 144 (<1%) 1,015 (4%) 6,509 (23%) 99.7 96.2 99.7 99.8 97.9 97.5 98.9 97.3 79.0 98.0 97.4 85.5 87.6 95.2 95.7 73.1 97.1 95.6 80.2 79.8 92.7 Baade et al. BMC Cancer (2015) 15:27 Page 4 of 13 Table 1 Cause-specific survival: invasive Melanoma in Queensland (1995–2008, follow-up to 2010) (Continued) Clarks level Into papillary dermis Filling papillary dermis Into reticular dermis Into subcutaneous fat Unknownb Positive nodes None At least one Metastasis No Yes 28,506 (99%) 148 (<1%) 99.3 68.7 95.3 39.4 92.9 27.6 28,421 (99%) 229 (<1%) 99.3 75.2 95.4 43.6 93.0 37.5 15,488 (54%) 6,164 (22%) 5,584 (19%) 715 (3%) 703 (2%) 99.9 99.6 98.0 93.4 92.4 99.4 95.3 86.6 70.0 80.9 98.5 92.7 80.3 60.7 73.9 Note: aModeled as a continuous variable in the prognostic model. b Actual values for missing data were estimated in the modeling process using multiple chained imputation. c sb = subungal; NOS = Not otherwise specified. assess and revise the adequacy of the final transformations of continuous covariates. Finally, we used a forward selection process to investigate whether the regression coefficients for any variables depended on follow-up time using the approach described by Royston & Parmar [20,21] with one degree of freedom for the time dependent covariate effects. This approach modifies the spline function of time used to model the baseline distribution function and is implemented using the stpm2t command in Stata. Covariates without time dependence were included in the model in standard fashion. Discrimination Development of the melanoma severity index (MSI) While a complete prognostic model containing all significant covariates may have its advantages, in particular being able to explain the greatest amount of variation in survival, a parsimonious model containing only important predictors may be preferred for clinical application [29]. To develop this reduced model, known as the MSI, we based the variable selection on the BIC statistics and the reduction in the per cent of explained variation when the variable was removed from the model [29,30]. For illustration, predicted survival probabilities and survival curves, along with their 95% confidence intervals were then calculated from the MSI based on specific combinations of the final set of prognostic variables. Validation and calibration The discrimination of a prognostic model reflects its ability to distinguish between patient outcomes, and is closely related to the idea of the proportion of variance that the model explains. We calculated Royston and Sauerbrei’s D statistic [28] as a measure of discrimination, and R2 as a measure of explained variation on the D natural scale of the model [28]. We also calculated the Harrel’s C discrimination index, which is a scored on a scale of 0 to 1. This can be taken to mean that if two cases are drawn at random, the c statistic is the probability that the person who survives the longest had the highest predicted survival. Values near 0.5 suggest the prognostic score is equivalent to a coin toss in determining which patient will live longer, while values near 0 or 1 indicate perfect discrimination. We assessed the importance of each variable in the prognostic model by examining the impact that removing the variable from the model had on D and R2 , both by removing just one of the variables, and D then by sequentially removing the variables from the model [28]. Calibration reflects prediction accuracy. A well-calibrated prognostic model assigns the correct mean survival probability at all levels of predicted risk. We used an internalexternal cross validation (IECV) method [31] for validating and assessing the calibration of the MSI. Briefly, the IECV method involves splitting Queensland into 9 geographical regions (Figure 1). The MSI is then fitted using data from eight of these regions. The linear predictor, Xβ was estimated and applied to both the fitted data (eight regions, “k”) and the excluded data (remaining 9th region, “(k)”). Values of Royston and Sauerbrei’s D statistic [28] were calculated from each result (Dk and D(k) respectively). If the predictive ability of the model is maintained, both values of D will be approximately equal. The difference between the two values (dk) was calculated, with its standard error being the square root of the sum of the squared standard errors for Dk and D(k). This process was repeated across the nine geographical regions. Baade et al. BMC Cancer (2015) 15:27 Page 5 of 13 Figure 1 Regions for internal-external validation, Queensland, Australia. Using these results, we assessed the calibration by comparing the predicted mean survival curves in each geographical region with the observed Kaplan-Meier survival estimates in that region. Results Patient characteristics Of 35,264 patients diagnosed with invasive melanoma in Queensland between 1995 and 2008, 34,384 were aged between 20 to 89 years at diagnosis, and 28,654 had no other histologically confirmed invasive melanoma diagnosed since 1982 (Table 1). The median follow up time, calculated as the median time to censoring, was 7.2 years. A total of 5,469 (19%) had died before the 31st December 2010. Almost a third of these deaths (31%) were ascribed to melanoma. Multivariable analysis the effect was small and of no clinical relevance, and so was not included in the final model. Therefore, the selected model had no time-dependent regression coefficients. Age at diagnosis was included in the model as a nontransformed continuous variable. Accurate modeling of the key predictor, thickness, was critical and was done in two stages. First, a smooth rank transform (SRT_thickness) was generated using the acd command in Stata [27] and further transformed for inclusion in the model (SRT_thickness2). Lack of fit was still apparent, so the original continuous thickness term was also included. The resulting analysis of martingale residuals showed satisfactory fit (results not shown). Concordance and discrimination All variables shown in Table 1 were considered in the initial prognostic model (see Figure 2). There was some evidence that sub-site had a time-dependent regression coefficient, so that the survival differential by subsite varied on the probit scale by follow up interval. However The effect of each variable on the discrimination (D) and explained variation ( R2 ) are shown in Table 2. D The final model explained 49% of the variation, with a D statistic of 1.55. Clearly, the majority of this performance was provided by the transformed thickness variables, with thickness alone accounting for nearly Baade et al. BMC Cancer (2015) 15:27 Page 6 of 13 Sex 1.0 0.9 0.8 1.0 0.9 0.8 Age at diagnosis 1.0 0.9 0.8 Morphology 1.0 0.9 0.8 Ulceration Survival probability Survival probability Survival probability Survival probability SSM LMM MM NOS 0 1 3 NM ALM MM in JN Other spec 0.7 0.6 0.5 0.4 0.3 0 1 Male 3 Female 5 7 10 Years since diagnosis 13 16 0.7 0.6 0.5 0.4 0.3 0 1 20-39 70-79 3 40-59 80-89 5 7 10 Years since diagnosis 13 60-69 0.7 0.6 0.5 0.4 0.3 0.7 0.6 0.5 0.4 0.3 No ulceration 0 1 3 Known ulceration 13 16 16 5 7 Years since diagnosis 13 16 5 7 Years since diagnosis Clarks Level 1.0 0.9 0.8 1.0 0.9 0.8 Thickness2 1.0 0.9 0.8 Subsite (a) 1.0 0.9 0.8 Subsite (b) Survival probability Survival probability Survival probability 0.7 0.6 0.5 0.4 Into papillary dermis Filling papillary dermis Into subc fat 0.7 0.6 0.5 0.4 <=0.50mm 1.01-1.50mm 2.01-4.00mm 0 1 3 0.51-1.00mm 1.51-2.00mm 4.01mm+ 13 16 0.7 0.6 0.5 0.4 0.3 0 1 Scalp chest/Axilla 3 Face/Lip/eyelid Abdomen/Hip 13 Ear/neck Survival probability 0.7 0.6 0.5 0.4 Back/buttocks Arm/Hand/fingers Shoulder Foot/heel/toes 0.3 0 1 Into reticular dermis 0.3 16 0.3 16 0 1 Leg/ankle 3 7 10 5 Years since diagnosis 13 5 7 Years since diagnosis 5 7 10 Years since diagnosis 3 5 7 10 Years since diagnosis 13 16 Positive nodes 1.0 0.9 0.8 1.0 0.9 0.8 Distant metastasis Survival probability Survival probability 0.7 0.6 0.5 0.4 0.3 0 1 None 3 Positive nodes 5 7 10 Years since diagnosis 13 16 0.7 0.6 0.5 0.4 0.3 0 1 No mets 3 Distant mets 13 16 5 7 10 Years since diagnosis Figure 2 Kaplan-Meier survival estimates by covariate group (1995–2008, follow up to 2010). Table 2 Effect of removing and adding variables from the prognostic model (Probit model with 3df) on discrimination (D) and explained variation (R2 ) D Removing variable Variable Clark’s level Morphology Gender Agec Positive lymph nodes Ulceration metastasis Body site Thicknessc a b Adding variable Cumulativeb Single variable only R2 D 0.476 0.473 0.469 0.458 0.451 0.438 0.420 0.401 D 1.108 0.643 0.277 0.398 0.990 0.998 1.068 0.303 1.307 R2 D 0.325 0.140 0.029 0.058 0.278 0.281 0.309 0.035 0.401 Harrell’s C 0.815 0.699 0.569 0.646 0.541 0.715 0.530 0.600 0.866 Added to thickness Dd 0.006 0.017 0.038 0.042 0.038 0.044 0.053 0.052 R2 d D 0.002 0.006 0.014 0.015 0.014 0.016 0.019 0.019 Order 1 2 3 4 5 6 7 8 9 D 1.521 1.513 1.500 1.467 1.457 1.409 1.358 1.307 - Singlea D 1.526 1.521 1.519 1.513 1.504 1.504 1.503 1.494 1.490 1.354 R2 D 0.478 0.476 0.475 0.473 0.470 0.470 0.470 0.467 0.466 0.419 Based on the prognostic regression model after removing only the given variable. Based on the prognostic regression model after sequentially removing the variables in the stated order (1 = first, 2 = second, ..). c Treated as transformed continuous variables (see text for details). The remaining variables are categorical. d Values represent the difference in fit statistics between the model with thickness only cand the model including thickness cand the shown variable. Baade et al. BMC Cancer (2015) 15:27 Page 7 of 13 41% of the variance and a D statistic of 1.32. However, when the remaining variables were combined, excluding thickness, they were able to explain similar amounts of the survival ( R2 ¼ 0:44; D ¼ 1:40) as thickness alone. D This result suggests relatively high correlation among the predictors. As variables entered as single variables in a model, the most discriminative variables were (in order) thickness, Clark’s level and metastasis. However, when removed from the full model (Table 2), Clark’s level (D: 0.005 and R2 : 0:002), morphology (D: 0.008 and D R2 : 0:003) and gender (D: 0.013 and R2 : 0:006) deD D creased the discrimination by small amounts and the explained variation by less than one percent. The BIC statistic was minimized by excluding Clark’s level and Morphology. Finally, Clark’s level, morphology and gender were removed to form the “parsimonious” model, hereby referred to as the Melanoma Severity Index (MSI). The Harrell-C statistic for the full model was 0.89 (95% CI = 0.890, 0.894), and 0.88 (0.877, 0.892) for the MSI. Generally the Harrell-C statistics for the individual variables were strongly correlated with the D and R2 statistics (Table 2), with the exception of metastaD ses and positive lymph nodes. Both these variables had low prevalence (Table 1). Additional unpublished sensitivity analyses suggested that the performance of the Harrell-C statistic in quantifying predictive ability of a variable depends on the prevalence of the feature in question. Final parameter estimates Internal-external validation The results of our internal-external cross validation (Figure 3, Table 4) suggests there is limited heterogeneity of the discrimination of the MSI model across the nine geographical regions of Queensland. The only exception was for Region 3, for which the model based on the rest of Queensland had higher discrimination when applied to this region. The mean survival curves predicted by the model generally agreed well with the observed (Kaplan-Meier) survival curves (Figure 3). Examples of model predictions for individual patients Examples of the estimated survival probabilities generated by the MSI for twelve hypothetical patients presenting with specific combinations of clinical characteristics and demographics are shown in Table 5 (and Figure 4), with comparison estimates generated by the AJCC Melanoma Patient Outcome Prediction Tool [32]. For example, the estimated 10-year survival of a 55 year old person diagnosed with a 0.45 mm thick melanoma on the back with no evidence of ulceration, positive lymph nodes or metastasis (Example 2) was 98.3%. In comparison, a the 10 year survival for a 75-year old person diagnosed with a 2.45 mm thick melanoma on the back with evidence of ulceration, but no positive lymph nodes or metastasis (Example 8) was 61.4%. When comparing the MSI and the MPOPT (Table 5), the predicted 10-year survival percentages for MSI were generally higher than for the MPOPT, with the exception of the two most advanced cases of melanoma. The parameter estimates from the full and the MSI multivariable probit regression model are shown in Table 3. The interpretation of the coefficients from this model is less familiar in a medical context than, for example, the coefficients from linear or logistic regression. A one-unit change in a covariate results in a one-beta change in risk on the probit (inverse normal probability) scale, where beta is the regression coefficient for the variable in question. However, in a more general sense, a positive beta coefficient means that an increase in the covariate raises the predicted probability of death from melanoma. Conversely, a negative beta coefficient means that an increase in the covariate reduces the predicted probability of death. The full model shows that the probability of death from melanoma was higher among males, older patients, those with thicker melanomas, lesions diagnosed on the scalp, abdomen/hip and foot, nodular melanomas, those with known ulceration, those with positive lymph nodes, evidence of distance metastasis and those with increasing level of invasion (Table 3). Similar direction and magnitude of effects for the included variables were observed for the MSI model (Table 3). Discussion This prognostic model for people diagnosed with invasive melanoma was developed using a large populationbased dataset containing information about the recognized prognostic features of melanoma with up to 16 years of follow-up. The flexible parametric survival model used for this development has many advantages over the standard Cox model used by previous studies, and the populationbased registry cohort removed the inherent biases associated with cohorts based on hospital clinic patients. In addition, by extracting additional clinical features from pathology forms, we were able to include clinical factors such as ulceration, number of positive lymph nodes and the presence of metastasis that are not typically available in these population-based data. There are several benefits of a prognostic model over and above those provided by standard Kaplan-Meier estimates. First, by incorporating demographic and clinical factors, the MSI output can be specifically targeted to an individual. Comparable Kaplan-Meier statistics would need to be generated from stratified subsamples, typically containing insufficient numbers to generate reliable estimates. Second, by incorporating the flexible parametric Baade et al. BMC Cancer (2015) 15:27 Page 8 of 13 Table 3 Multivariable parameter estimates from the Full and MSI modelsa Variable Gender Male Female Age at diagnosis (transformed) Thickness (original) SRT_Thickness (transformed)b Body site Scalp Face/Lip/Eyelid Ear/Neck Chest/Axilla Abdomen + Hip Back/Buttocks Upper arm/Forearm/Hand/Finger/sb.hand Foot/sb.foot/heel/toec Thigh/leg/ankle Shoulder Morphology SSM Nodular melanoma MM in junctional naevus Lentigo Maligna melanoma Acral lentiginous melanoma Other specified melanoma Malignant melanoma NOSd Ulceration No ulceration Known ulceration Positive lymph nodes None At least one Metastasis No Yes Clarks Level Into papillary dermis Filling papillary dermis Into reticular dermis Into subcutaneous fat a Beta coefficient [95% confidence interval] Full model 0.000a −0.147 [−0.207, −0.086] 0.007 [0.005, 0.008] 0.044 [0.030, 0.058] 1.593 [1.394, 1.792] 0.007 [0.005, 0.009] 0.038 [0.025, 0.050] 1.830 [1.697, 1.963] MSI model Not included 0.376 [0.233, 0.519] −0.063 [−0.173, 0.046] −0.025 [−0.141, 0.091] −0.108 [−0.250, 0.033] 0.148 [−0.010, 0.307] 0.000a c 0.348 [0.210, 0.486] −0.110 [−0.215, −0.004] −0.041 [−0.156, 0.073] −0.102 [−0.244, 0.039] 0.143 [−0.015, 0.301] 0.000a −0.361 [−0.463, −0.259] 0.089 [−0.089, 0.267] −0.291 [−0.376, −0.206] −0.099 [−0.207, 0.009] −0.321 [−0.425, −0.218] 0.174 [−0.031, 0.378] −0.251 [−0.339, −0.163] −0.079 [−0.187, 0.029] 0.000a 0.076 [−0.005, 0.158] 0.085 [−0.154, 0.325] −0.084 [−0.244, 0.077] −0.081 [−0.402, 0.241] −0.232 [−0.355, −0.109] 0.014 [−0.057, 0.086] 0.000a 0.357 [0.277, 0.437] 0.000a 0.753 [0.593, 0.912] 0.000a 1.062 [0.861, 1.262] 0.000a 0.153 [0.052, 0.253] 0.193 [0.072, 0.314] 0.197 [0.024, 0.370] Not included 0.000a 0.378 [0.300, 0.456] 0.000a 0.779 [0.620, 0.938] 0.000a 1.062 [0.861, 1.262] Not included Note: denotes the reference category. b SRT_Thickness is the smooth rank transform of the thickness variable (see Methods for details). c sb = subungal; dNOS = Not otherwise specified. Baade et al. BMC Cancer (2015) 15:27 Page 9 of 13 Region 1 1.00 Region 2 A 1.00 Region 3 A 1.00 A B C Survival probability Survival probability 0.90 0.80 0.70 0.60 B C 0.90 0.80 0.70 0.60 B C Survival probability 0.90 0.80 0.70 0.60 D 0.50 0.40 0.30 0 2 4 6 8 10 D 0.50 0.40 0.30 0 2 4 6 8 10 D 0.50 0.40 0.30 0 2 4 6 8 10 Years from diagnosis Years from diagnosis Years from diagnosis Region 4 1.00 Region 5 A 1.00 Region 6 A 1.00 A B C Survival probability Survival probability 0.90 0.80 0.70 0.60 B C 0.90 0.80 0.70 0.60 B C Survival probability 0.90 0.80 0.70 0.60 D 0.50 0.40 0.30 0 2 4 6 8 10 D 0.50 0.40 0.30 0 2 4 6 8 10 D 0.50 0.40 0.30 0 2 4 6 8 10 Years from diagnosis Years from diagnosis Years from diagnosis Region 7 1.00 Region 8 A 1.00 Region 9 A 1.00 A B C Survival probability Survival probability 0.90 0.80 0.70 0.60 B C 0.90 0.80 0.70 0.60 B C Survival probability 0.90 0.80 0.70 0.60 D 0.50 0.40 0.30 0 2 D 0.50 0.40 0.30 D 0.50 0.40 0.30 Years from diagnosis 4 6 8 10 0 2 Years from diagnosis 4 6 8 10 0 2 Years from diagnosis 4 6 8 10 Observed Kaplan-Meier survival in excluded region Predicted survival in excluded region Figure 3 Predicted and observed survival curves by region from the internal–external cross-validation approach. A: < 70th centile B: 70-80th centile C: 80-90th centile D: >90th centile approach to the prognostic model, it provides the ability to generate other estimates of survival, including longer-term survival, conditional survival and the predicted life expectancy following a melanoma diagnosis. The single most important feature in the prognostic model was thickness. This has been shown to be an important predictor in previous models, [11,13,32,33] even when restricted to localized [13] or thin [2] melanomas. The MSI can be used to highlight the importance of detecting melanoma early from a public health perspective by quantifying the expected reduction in survival as the thickness at diagnosis increases. While this prognostic model has not yet been validated against an external, independent dataset, we Table 4 Evaluation of heterogeneity of the measure of separation (Discrimination, D) across geographically defined regions: internal–external cross-validation using the MSI model Region (k)1 1 2 3 4 5 6 7 8 9 1 Nk 1,294 1,438 919 1,574 1,461 2,313 3.495 10,578 3,982 Melanoma deathsk 84 80 54 91 104 138 216 601 212 Other deathsk 150 151 101 192 204 287 446 1,376 489 D(k)Omitting region k 1.50 1.50 1.49 1.49 1.50 1.49 1.53 1.52 1.53 DkPredicted in region k 1.57 1.51 2.34 1.73 1.52 1.75 1.38 1.46 1.40 dk(= Dk– D(k)) 0.07 0.01 0.85 0.24 0.01 0.26 −0.15 −0.06 −0.13 SE (dk) 0.14 0.14 0.24 0.14 0.13 0.12 0.08 0.06 0.09 Cases for which the geographical region was undefined (n = 1,600, 120 melanoma deaths, 373 other deaths) were retained in the validation process. Baade et al. BMC Cancer (2015) 15:27 Page 10 of 13 Table 5 Predicted 10-year survival percentages for twelve hypothetical melanoma patients using the MSI, including comparisons with the AJCC Melanoma Patient Outcome Prediction Tool (MPOPT)a Characteristic Age 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 35 55 55 75 55 75 75 75 75 85 75 75 Thick (mm) 0.45 0.45 1.00 1.00 1.45 1.45 2.45 2.45 4.45 4.45 4.45 4.45 Site Back Back Back Back Back Back Back Back Scalp Scalp Scalp Scalp Ulceration No No No No No No No Yes Yes Yes Yes Yes Positive lymph nodes No No No No No No No No No No Yes Yes Mets No No No No No No No No No No No Yes Estimated survival (%) after ten years MSI 98.8 [99–99] 98.3 [98–99] 91.7 [91–93] 89.4 [88–91] 86.2 [85–88] 83.0 [81–85] 74.8 [72–77] 61.4 [58–64] 38.9 [34–44] 36.3 [31-42] 14.5 [10–20] 1.7 [1–3] MPOPTa 97.5 [97–99] 97.5 [97–99] 90.0 [88–92] 82.8 [78–88] 83.7 [82–86] 69.6 [63–77] 59.6 [51–69] 42.0 [33–54] 35.2 [27–47] 35.2 [27–47] 29.1 [22–39] 15.9 [6–44] Notes: aFor comparison with the MPOPT (www.melanomaprognosis.org), Scalp and Back were categories as “Axial”, and positive lymph nodes or metastasis were both considered indicative of “Regional metastasis). 1 Survival (%) 20 40 60 80 100 Survival (%) 20 40 60 80 100 2 Survival (%) 20 40 60 80 100 3 Survival (%) 20 40 60 80 100 4 0 0 0 0 0 2 4 6 8 Years since diagnosis 10 0 6 8 2 4 Years since diagnosis 10 0 2 4 6 8 Years since diagnosis 10 0 6 8 2 4 Years since diagnosis 10 5 Survival (%) 20 40 60 80 100 Survival (%) 20 40 60 80 100 6 Survival (%) 20 40 60 80 100 7 Survival (%) 20 40 60 80 100 8 0 0 0 0 0 0 2 4 6 8 Years since diagnosis 10 0 2 4 6 8 Years since diagnosis 10 0 2 4 6 8 Years since diagnosis 10 2 4 6 8 Years since diagnosis 10 9 Survival (%) 20 40 60 80 100 Survival (%) 20 40 60 80 100 10 Survival (%) 20 40 60 80 100 11 Survival (%) 20 40 60 80 100 12 0 0 0 0 0 2 4 6 8 Years since diagnosis 10 0 8 2 4 6 Years since diagnosis 10 0 8 2 4 6 Years since diagnosis 10 0 6 8 2 4 Years since diagnosis 10 Figure 4 Predicted 10-year mean survival curves, for twelve hypothetical melanoma patients. Baade et al. BMC Cancer (2015) 15:27 Page 11 of 13 have examined the internal reliability by assessing its consistency across a variety of different geographical areas within the state, areas that are characterized by differing health patterns and life expectancy. These unmeasured factors are not considered by the prognostic model, so the consistency of the model’s discriminatory performance across these regions is encouraging. Validation using an external dataset would entail calculating the prognostic index using the parameter estimates from this study cohort applied to the covariate values of the secondary dataset. Similarly, the baseline survival function is calculated in the external dataset using both the parameter vector from this study cohort, together with the (log) time values in the external dataset and the set of spline knots used in the current cohort [34]. We did not have information on mitotic rate, and thus have not been able to examine the effect this had on melanoma survival. Mitotic rate is recognized as an important prognostic factor in the final version of the 2009 AJCC Melanoma Staging and Classification [11]. That this variable was omitted here is a limitation. In addition there are other factors that have been shown to influence melanoma survival that have not been included in this model [33], such as HIV infection, race and socioeconomic status. We also did not have information about concomitant diseases that may have independently lowered predicted survival. Importantly, we did include variables that were demonstrated by others to be significant predictors of melanoma survival by other studies, notably thickness, age, body site and ulceration [13,32]. The use of multiple imputation for missing data assumed that the data are missing at random (MAR). We were not able to rule out that the data are missing not at random (MNAR), and it remains possible that there is some unmeasured characteristic of the treating clinician or pathology laboratory that impacted on the completeness of the registry data. Finally, since we included only patients with one diagnosed melanoma to ensure a greater link between the melanoma characteristics and survival outcome, the predicted survival outcomes would no longer be relevant if a subsequent melanoma was diagnosed. Given the variation in cohort selection, clinical characteristics and statistical methods, direct comparisons with the results of our study to those published recently [2,13,32] for melanoma are difficult. In particular, the comparisons with the AJCC Melanoma Patient Outcome Prediction Tool (United States) [13,32] demonstrate there are important differences between the two countries. It is unknown whether this is due to the statistical method or to cohort selection, since the MPOPT is based on patients selected from major cancer centres and clinical trial cooperative groups. Alternatively, differences in diagnostic or management practices may have led to important differences in the survival outcomes expected by people diagnosed with melanoma in the two countries. From a clinical perspective, this MSI could be readily applied to individual melanoma patients. We plan to incorporate a MSI online dissemination tool into a broader package for primary care physicians to assist in discussions about prognosis in the clinical setting. One of the many advantages of using the flexible parametric modeling approach is that these models can be readily extended to consider estimates of conditional survival [20] and loss of life expectancy [35] due to their diagnosis in comparison to the general population. We have previously shown that people diagnosed with localized or regional melanoma have negligible excess mortality compared to the general population once they have survived ten years post-diagnosis [36,37]. Such models enable more precise estimates of conditional survival by incorporating clinically relevant factors in the estimation. Further work is also planned on the impact of competing risks within the flexible parametric framework [20,38]. We did not find any time-dependent regression coefficients of sufficient impact to include them in the final model. Proponents of other prognostic models using large epidemiological cancer datasets have observed time-varying coefficients [38-40] on a hazard scale. However, in our probit model, hazard ratios comparing any two values of a covariate can be shown to tend toward 1 as follow-up time increases. Thus, all hazard ratios are time-dependent. When there are multiple timedependent coefficients, interpreting the time dependent hazard ratios can be difficult in the log cumulative hazard framework of the Royston-Parmar models [20]. The reason is that hazard ratios depend on the values of more than just one covariate. Alternatives, including modelling on the log excess hazard scale, may offer more interpretable options when time dependent coefficients are present [41]. Conclusions The MSI serves to identify and weight key parameters of prognostic importance in melanoma, and therefore allow a finer gradation of the prognostic value of interventions than the simple dichotomous variable of death due to melanoma, with greater statistical power for comparisons. It also enables predictions of survival outcomes soon after the intervention, rather having a 5 or ten year delay for follow-up. Thus, beyond the immediate clinical use, the MSI may have important public health and research applications for evaluations of public health interventions aimed at reducing deaths from melanoma. Abbreviation MSI: Melanoma severity index. Baade et al. BMC Cancer (2015) 15:27 Page 12 of 13 Competing interests The authors declare that they have no competing interests. 13. Authors’ contributions JA, AG, MW, PB, PY conceived the study, PB and PR carried out the statistical analysis, all authors contributed to the drafting and revision of the manuscript, all authors read and approved the final manuscript. Acknowledgement We are indebted to the staff at the Queensland Cancer Registry for collecting and maintaining these data. Financial support Prof Baade was supported by an NHMRC Career Development Fellowship (ID1005334). Dr Philippa Youl was supported by an NHMRC Early Career Fellowship (ID1054038). The funder had no input into the design and conduct of the study, collection, management, analysis, and interpretation of the data, preparation, review, or approval of the manuscript or the decision to submit the manuscript for publication. Author details 1 Cancer Council Queensland, 553 Gregory Terrace, Fortitude Valley, Queensland, Australia. 2School of Public Health and Social Work, Queensland University of Technology, Brisbane, Queensland, Australia. 3Griffith Health Institute, Griffith University, Gold Coast, Queensland, Australia. 4MRC Clinical Trials Unit at University College London, London, UK. 5Dermatoepidemiology Unit, V A Medical Center, Providence, RI, USA. 6Department of Dermatology, Rhode Island Hospital, Providence, RI, USA. 7Departments of Dermatology and Epidemiology, Brown University, Providence, RI, USA. 8Harvard School of Public Health, Harvard University, Boston, MA, USA. Received: 18 March 2014 Accepted: 14 January 2015 22. References 1. GLOBOCAN 2008, Cancer incidence and mortality worldwide: IARC CancerBase No. 10 [Internet] [http://globocan.iarc.fr/] 2. Green AC, Baade P, Coory M, Aitken JF, Smithers M. Population-based 20-year survival among people diagnosed with thin melanomas in Queensland, Australia. J Clin Oncol. 2012;30(13):1462–7. 3. Hu CY, Xing Y, Cormier JN, Chang GJ. Assessing the utility of cancerregistry-processed cause of death in calculating cancer-specific survival. Cancer. 2013;119(10):1900–7. 4. Eisemann N, Jansen L, Holleczek B, Waldmann A, Luttmann S, Emrich K, et al. Up-to-date results on survival of patients with melanoma in Germany. Br J Dermatol. 2012;167(3):606–12. 5. Jeffreys M, Sarfati D, Stevanovic V, Tobias M, Lewis C, Pearce N, et al. Socioeconomic inequalities in cancer survival in New Zealand: the role of extent of disease at diagnosis. Cancer Epidemiol Biomarkers Prev. 2009;18(3):915–21. 6. AIHW. Cancer survival and prevalence in Australia: period estimates from 1982 to 2010. Asia Pac J Clin Oncol. 2013;9(1):29–39. 7. Forsea AM, Del Marmol V, de Vries E, Bailey EE, Geller AC. Melanoma incidence and mortality in Europe: new estimates, persistent disparities. Br J Dermatol. 2012;167(5):1124–30. 8. Gimotty PA, Guerry D, Ming ME, Elenitsas R, Xu X, Czerniecki B, et al. Thin primary cutaneous malignant melanoma: a prognostic tree for 10-year metastasis is more accurate than American joint committee on cancer staging. J Clin Oncol. 2004;22(18):3668–76. 9. Xing Y, Chang GJ, Hu CY, Askew RL, Ross MI, Gershenwald JE, et al. Conditional survival estimates improve over time for patients with advanced melanoma: results from a population-based analysis. Cancer. 2010;116(9):2234–41. 10. Steyerberg EW, Moons KG, van der Windt DA, Hayden JA, Perel P, Schroter S, et al. Prognosis research strategy (PROGRESS) 3: prognostic model research. PLoS Med. 2013;10(2):e1001381. 11. Balch CM, Gershenwald JE, Soong S-j, Thompson JF, Atkins MB, Byrd DR, et al. Final version of 2009 AJCC melanoma staging and classification. J Clin Oncol. 2009;27(36):6199–206. 12. Balch CM, Soong S-J, Gershenwald JE, Thompson JF, Reintgen DS, Cascinelli N, et al. Prognostic factors analysis of 17,600 melanoma patients: validation 23. 24. 25. 14. 15. 16. 17. 18. 19. 20. 21. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. of the American joint committee on cancer melanoma staging system. J Clin Oncol. 2001;19(16):3622–34. Soong SJ, Ding S, Coit D, Balch CM, Gershenwald JE, Thompson JF, et al. Predicting survival outcome of localized melanoma: an electronic prediction tool based on the AJCC Melanoma Database. Ann Surg Oncol. 2010;17(8):2006–14. Egger ME, Dunki-Jacobs EM, Callender GG, Quillo AR, Scoggins CR, Martin 2nd RC, et al. Outcomes and prognostic factors in nodular melanomas. Surgery. 2012;152(4):652–9. discussion 659–660. Coory M, Smithers M, Aitken J, Baade P, Ring I. Urban–rural differences in survival from cutaneous melanoma in Queensland. Aust N Z J Public Health. 2006;30(1):71–4. Lyth J, Hansson J, Ingvar C, Mansson-Brahme E, Naredi P, Stierner U, et al. Prognostic subclassifications of T1 cutaneous melanomas based on ulceration, tumour thickness and Clark’s level of invasion: results of a population-based study from the Swedish Melanoma Register. Br J Dermatol. 2013;168(4):779–86. Lange JR, Palis BE, Chang DC, Soong SJ, Balch CM. Melanoma in children and teenagers: an analysis of patients from the National Cancer Data Base. J Clin Oncol. 2007;25(11):1363–8. Grulich AE, Swerdlow AJ, Silva IDS, Beral V. Is the apparent rise in cancer mortality in the elderly real? analysis of changes in certification and coding of cause of death in England and Wales, 1970–1990. Int J Cancer. 1995;63(2):164–8. Schemper M, Smith TL. A note on quantifying follow-up in studies of failure time. Control Clin Trials. 1996;17(4):343–6. Royston P, Lambert PC. Flexible parametric survival analysis using stata: beyond the cox model. College Station, Texas: Stata Press; 2011. Royston P, Parmar MK. Flexible parametric proportional-hazards and proportional-odds models for censored survival data, with application to prognostic modelling and estimation of treatment effects. Stat Med. 2002;21(15):2175–97. Rubin D. Multiple imputation for nonresponse in surveys. New York: Wiley; 1987. White I, Royston P, Wood A. Multiple imputation using chained equations: issues and guidance for practice. Stat Med. 2011;30(4):377–99. Reid N. A conversation with Sir David Cox. Stat Sci. 1994;9(3):439–55. Sauerbrei W, Royston P. Building multivariable prognostic and diagnostic models: transformation of the predictors by using fractional polynomials. J R Stat Soc (Series A). 1999;162:71–94. Royston P, Sauerbrei W. Building multivariable regression models with continuous covariates in clinical epidemiology–with an emphasis on fractional polynomials. Methods Inf Med. 2005;44(4):561–71. Royston P: A smooth covariate rank transformation for use in regression models with a sigmoid dose–response function. Submitted. Royston P, Sauerbrei W. A new measure of prognostic separation in survival data. Stat Med. 2004;23(5):723–48. Ambler G, Brady AR, Royston P. Simplifying a prognostic model: a simulation study based on clinical data. Stat Med. 2002;21(24):3803–22. Royston P, Sauerbrei W. Multivariable model - building : a pragmatic approach to regression anaylsis based on fractional polynomials for modelling continuous variables. Hoboken: Wiley; 2008. Royston P, Parmar MK, Sylvester R. Construction and validation of a prognostic model across several studies, with an application in superficial bladder cancer. Stat Med. 2004;23(6):907–26. Balch CM, Soong SJ, Gershenwald JE, Thompson JF, Coit DG, Atkins MB, et al. Age as a prognostic factor in patients with localized melanoma and regional metastases. Ann Surg Oncol. 2013;20(12):3961–8. Wisco OJ, Sober AJ. Prognostic factors for melanoma. Dermatol Clin. 2012;30(3):469–85. Royston P, Parmar M, Altman DG: External validation and updating of a prognostic survival model (Research report No. 307). In.: Department of statistical science, University College London; 2010. Andersson TM-L, Dickman P, Eloranta S, Lambe M, Lambert P: Estimating the loss in expectation of life due to cancer using flexible parametric survival models. Stat Med; In press. Baade PD, Youlden DR, Chambers SK. When do I know I am cured? Using conditional estimates to provide better information about cancer survival prospects. Med J Aust. 2011;194(2):73–7. Yu XQ, Baade PD, O'Connell DL. Conditional survival of cancer patients: an Australian perspective. BMC Cancer. 2012;12:460. Baade et al. BMC Cancer (2015) 15:27 Page 13 of 13 38. Hinchliffe SR, Lambert PC. Flexible parametric modelling of cause-specific hazards to estimate cumulative incidence functions. BMC Med Res Methodol. 2013;13:13. 39. Cramb SM, Garvey G, Valery PC, Williamson JD, Baade PD. The first year counts: cancer survival among Indigenous and non-Indigenous Queenslanders, 1997–2006. Med J Aust. 2012;196(4):270–4. 40. Colzani E, Liljegren A, Johansson AL, Adolfsson J, Hellborg H, Hall PF, et al. Prognosis of patients with breast cancer: causes of death and effects of time since diagnosis, age, and tumor characteristics. J Clin Oncol. 2011;29(30):4014–21. 41. Crowther MJ, Lambert PC. A general framework for parametric survival analysis. Stat Med. 2014;33(30):5280–97. Submit your next manuscript to BioMed Central and take full advantage of: • Convenient online submission • Thorough peer review • No space constraints or color figure charges • Immediate publication on acceptance • Inclusion in PubMed, CAS, Scopus and Google Scholar • Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit