1. Background
Development of risk prediction models is currently a common and appealing area of research (1-3). The discovery and study of new biomarkers has been a driving force in clinical medicine in the last two decades. The development of prediction algorithms based on multivariate regression models began to show promise several decades ago. Predictive models are attractive in that they allow individuals to use their risk factors profile for calculating their corresponding risk of developing the event of interest in the future (3).
Parallel with the development of predictive models, biomarker research started to emerge on an impressively large scale (1). In light of the pace at which new risk markers are being discovered by scientists, statisticians and clinicians face challenges and opportunities regarding how best to evaluate these biomarkers and how to incorporate them into new prediction models. Furthermore, many predictive models yield risk values that fall into the intermediate range for persons (3). While decisions are easiest when the estimated risk is very low or very high, enhancement to the extant models has been sought to increase their predictive performance.
Toward this aim, new biomarkers have been added to relevant models (4). Recently, researchers have asserted that the predictive performance of one model is superior to another, an assertion that frequently has been challenged by statistical reviewers of scientific journals asking for rigorous statistical justification for researchers’ statements. Demonstration of a statistically significant association of a new biomarker with cardiovascular risk is not sufficient (5-8). What is the best method for quantifying the improvement in risk prediction offered by these new models? The answer to this question plays a pivotal role in dealing with new risk markers. This question, of course, conceives with the more basic concept of how to assess the performance of a risk prediction model (9, 10). Researchers take advantage of a host of methods and metrics in order to assess the performance of prediction models. In longitudinal survival analysis, a concordance measure statistic has been developed based on the same concept as the receiver operating characteristic curve. Several new measures have recently been proposed:
1. Reclassification tables, net reclassification improvement (NRI), and integrated discrimination improvement (IDI) for binary outcome and survival analysis (1, 2, 6-9, 11).
2. A variant of the Hosmer-Lemeshow X2 (Nam-D’Agostino X2) for measuring the calibration of a survival model (12).
3. A variant of the AROC (Harrell’s C statistic) for measuring the discrimination of a survival model (12-15).
Commonly-used, user-friendly statistical packages such as SPSS have not provided calculations for novel predictive performances statistics. Many studies thus either take no notice of the novel statistics or inadvertently use novel analyses, i.e., using original NRI while dealing with censored survival data. Although open-source, free statistical packages can be used to calculate novel statistics, they need some knowledge of programming, which renders their usage limited. Furthermore, predictive models in medicine are usually developed based on the regression model, which is simply a line of command in Stata, while requiring hundreds of lines in R. There are currently no commands available in Stata by which to calculate Nam-D’Agostino X2, IDI, or NRI. Therefore, we wrote the adpredsurv command for calculating IDI and NRI for survival models and the ndacs command for calculating Nam-D’Agostino X2 for survival models.
2. Objectives
User-friendly statistical software capable of implementing novel statistical procedures is conspicuously lacking. This shortcoming has restricted implementation of novel model assessment methods. We aimed to provide Stata commands to help researchers obtain the aforementioned statistical indices.
3. Materials and Methods
3.1. Design
3.1.1. Added Predictive Ability for Survival Models
3.1.1.1. IDI
IDI has previously been defined by Polak et al. as for logistic models (16). IDI for survival-based regression models is defined as follows.
Equation 1.
This ratio can be interpreted as the proportion of variance explained by the model. Chambless et al. (17) have labeled this ratio R2(t) and have shown that 0 ≤ R2(t) ≤ 1. The parameter can be estimated for fitted survival Equation:
Equation 2.
The IDI thus can be interpreted as the difference in the variance explained by different models:
Equation 3.
3.1.1.2. NRI
Pencina et al. have reformulated NRI so that it can be interpreted prospectively (18):
Equation 4.
After sum simplification, the above will reduce to:
Equation 5.
The extension of this formula to survival analysis is immediate, with P (event), P (event, up), P (event, down) all estimated using the Kaplan-Meier approach (18).
Assuming that out of a total of n individuals, nU are reclassified upward and nD downward, can be written as:
Equation 6.
The Formula 4 does not depend on the number and even the existence of categories as. In fact, in this way, the NRI is calculated based on the presumption that probabilities of the event of interest among those reclassified downward or upward is obtained by pooling all individuals who share the same reclassification (18).
3.1.1.3. Nam-D’Agostino X2
If the predicted risk from a survival model is divided into j groups, j = 1, 2, …, M, KMj is the Kaplan-Meier estimated observed incidence rate for the event of interest in jth group, and is the average model-based predicted probability of the event of interest, then the Nam-D’Agostino X2 is estimated as:
Equation 7.
3.1.2. Commands
3.1.2.1. “Adpredsurv” Command
3.1.2.1.1. Syntax
Adpredsurv outcome old risk new risk, cut point (numlist) end time (numlist)
3.1.2.1.2. Description
“Adpredsurv” calculated relative and absolute IDI as well as cutpoint-based and cutpoint-free NRI for survival models.
3.1.2.1.3. Options
“Cut point (numlist)” specifies the numbers that present cut points of risk based on the old model on which the new model is to be evaluated.
“End time (numlist)” specifies the survival time at which the analysis is to be performed.
3.1.2.2. “Ndacs” Command
3.1.2.2.1. Syntax
Ndacsriskvar timevar, group(numlist) end(numlist).
“Ndacs” calculates Nam-D’Agostino X2 (and is P value) which is a version of the Hosmer-Lemeshow test for calibration and Harrell’s C index of discrimination (and its 95% confidence intervals) for survival models. The risk estimated based on the survival model under the investigation is specified by “riskvar”. Survival time is specified in “timevar”.
3.1.2.2.2. Options
“group()” specifies the number of quantiles to be used to group the data for the Nam-D’Agostino goodness of fit test. “group (10)” is typically specified. “end()” specifies the survival time at which the analysis is to be performed. It is typically set at the median follow-up time.
3.2. Example
3.2.1. Study Population
We used a real data set of the Tehran Lipid and Glucose Study (TLGS) to predict incidence of cardiovascular diseases (CVD). Detailed descriptions of the TLGS have been reported elsewhere (19). In brief, the TLGS is a large-scale, long-term, community-based prospective study performed on a representative sample of residents of District 13 of Tehran, the capital of Iran. Age and sex distributions of the population in the district were representative of the overall population of Tehran at the time of the baseline examination. The TLGS has two major components: a cross-sectional prevalence study of non-communicable diseases and associated risk factors implemented between March 1999 and December 2001, and a prospective follow-up study. Data collection is ongoing, designed to continue for at least 20 years, at three-year intervals. A total of 27,340 residents aged ≥ 3 years were invited by telephone call, of which 15,010 residents participated in the first examination cycle and another 3,551 residents were first examined in the second examination cycle. Participants were categorized into the cohort (n = 10 394) and intervention groups (n = 8 167), the latter to be educated for implementation of lifestyle modifications. For the current study, among participants aged ≥30 years (n = 9 752), we selected those who participated in the follow-up study until March 20, 2009 (n = 8 795). We used data on 4,052 women with complete data on covariates, contributing to a 42 659 person-year follow up. At the time of this study, the median follow-up time was 11.5 years.
3.2.2. Clinical and Laboratory Measurements
Information was collected by a trained interviewer using a pretested questionnaire and included demographic data, drug history, past medical history of CVD, hypertension, and diabetes and smoking status (20). After a 15-minute rest in a seated position, using a standardized mercury sphygmomanometer (calibrated by the Iranian institute of standards and industrial researches), two measurements of blood pressure were taken on the right arm. The mean of the two measurements was considered the participant’s blood pressure.
After fasting overnight for 12 to 14 hours, a blood sample was drawn between 7:00 a.m. and 9:00 a.m. from all study participants. All samples were analyzed when internal quality control met the acceptable criteria (19). All the blood analyses were undertaken on the day of blood collection at the TLGS research laboratory. Total cholesterol (TC) was assayed, using the enzymatic colorimetric method with cholesterol esterase and cholesterol oxidase. High-density lipoprotein cholesterol (HDL-C) was measured after precipitation of the apolipoprotein B containing lipoproteins with phosphotungistic acid. Analyses were performed using Pars Azmon kits (Pars Azmon Inc., Tehran, Iran) and a Selectra 2 auto-analyzer (Vital Scientific, Spankeren, Netherlands). An enzymatic colorimetric method with glucose oxidase was used for the measurement of plasma glucose. The standard two-hour post-challenge plasma glucose (2 hours-PCPG) test was used for participants not on glucose-lowering agents.
4. Results
4.1. Outcome Measurements
Details of cardiovascular outcomes have been published elsewhere (21). In the follow-up study, telephone calls are made by a trained nurse to every TLGS participant. Participants are asked about any medical event during the previous year or whether related events have occurred. When appropriate, a trained physician collects complementary data during a home visit and or a visit to the respective hospital to collect data from the participants’ medical records. In the case of mortality, data are collected from the hospital or the death certificate by an authorized local physician. Collected data are evaluated by an outcome committee consisting of a principal investigator, an internist, an endocrinologist, a cardiologist, an epidemiologist, and the physician who collects the outcome data. Other experts are invited for evaluation of non-communicable disorders on an as-needed basis.
A specific outcome for each event is assigned according to international statistical classification of diseases and related health problems CRITERIA, 10th Revision, and the American heart association classification for cardiovascular events (19, 22, 23). Coronary heart disease (CHD) includes cases of definite myocardial infarction (MI) diagnosed by electrocardiogram (ECG) and biomarkers, probable MI (positive ECG findings plus cardiac symptoms or signs and biomarkers showing negative or equivocal results), unstable angina pectoris (new cardiac symptoms or changing symptom patterns and positive ECG findings with normal biomarkers), angiographic-proven CHD, and CHD death. CVD is specified as a composite measure of any CHD events, stroke, or cerebrovascular death.
4.2. Definition of Terms
Current smoker was defined as a person who smokes cigarettes daily or occasionally. A previous history of CVD reflected any prior diagnosis of CVD by a physician. In accordance with the definition provided by the American diabetes association, participants were classified as having diabetes at the baseline if they met at least one of these criteria: FPG ≥ 7 mmolL-1, or 2 hours-PCPG ≥ 11.1 mmolL-1, or taking anti-diabetic medication (19). The family history of premature CVD was obtained by asking participants whether any member in their immediate family (first-degree relatives) had experienced a fatal or nonfatal MI, stroke, or sudden cardiac arrest. The event was considered premature if it occurred before the age of 55 years in male relatives and before 65 in female relatives.
4.3. Statistics Analysis
We used the Weibull regression model for analysis of outcomes (CVD). Baseline Weibull regression model was developed based on the traditional risk markers (e.g., age, smoking, systolic blood pressure, use of anti-hypertensive drugs, total and HDL cholesterol, and diabetes). An improved Weibull regression model was developed by adding family history of premature CVD to the baseline model to the basic Weibull model.
We set the statistical significance level at a two-tailed type I error of 0.05 and used Stata version 12.0 (StataCorp, College Station, Texas USA) for all statistical analyses.
It is certified that all applicable governmental and institutional regulations regarding the ethical use of human volunteers were followed during this research. The ethical committee of the Research Institute for Endocrine Sciences approved this study. Informed written consent was obtained from all participants. The investigations reported herein have been carried out in accordance with the principles of the declaration of Helsinki as revised in 2008.
4.4. Assessment of Model Performance
We used several criteria to compare the overall predictive values of alternative models.
5. Discussion
In the survival analysis, discrimination capacity is quantified by Harrell’s C statistic. It calculates the probability that for each randomly selected tie of two persons who did and did not develop the event of interest, a person who developed an event of interest at a certain specific time has a higher risk score than a randomly selected person who did not develop an event during the same specific follow-up interval (15). The 95% confidence intervals for Harrell’s C statistic of different models were estimated with bootstrap resampling.
5.1. Calibration
Calibration quantified the magnitude to which predicted probabilities closely agree numerically with actual outcomes. Calibration was examined by calculating the Nam-D’Agostino X2 and using a test similar to the Hosmer-Lemeshow test as suggested in the literature (12).
5.1.1. Added Predictive Capacity
Absolute and relative IDI and cut-point-based and cut-point-free NRI were used as measures of predictive ability added to the baseline survival-based regression model by paraclinical parameters (2). The bootstrapping method was implemented in order to obtain 95% confidence intervals (95% CIs). Tables 1 and 2 presents the novel analysis obtained from Stata.
Table 1. Novel Analysis
| Basic Model | Enhanced Model | |
|---|---|---|
| Nam-D’Agostino X2 | 57.14 | 41.2 |
| P value | < 0.001 | < 0.001 |
| Harrell’s C index (95% CIs) | 0.832 (0.812 - 0.852) | 0.837 (0.817 - 0.856) |
Table 2. Predictive Performances of the Basic Framingham’s “general CVD risk” Algorithm vs. Enhanced Modela
| Added Predictive Values | 95% CIs | P Value |
|---|---|---|
| Absolute IDI | 0.0047 (-0.0004 - 0.0099) | 0.073 |
| Relative IDI | 0.0506 (-0.0056 - 0.1067) | 0.078 |
| Cutpoint-based NRI | 0.0142 (-0.0285 - 0.0569) | 0.514 |
| Cutpoint-free NRI | 0.1697 (0.0608 - 0.2785) | 0.002 |
aIDI, integrated discriminatory improvement index; NRI, net reclassification improvement index (95% CIs).
The Framingham’s “general CVD risk” algorithm incorporated age, systolic blood pressure, use of blood pressure lowering drugs, total and high-density lipoprotein cholesterol, smoking, and diabetes. The enhanced model was developed by adding a family history of premature CVD to the basic Framingham’s “general CVD risk” algorithm.
5.1.2. Conclusion
In all major fields of modern medicine (e.g., cardiovascular disease, cancer, and diabetes) risk prediction models continue to emerge (24-34). The predictive performance of these models needs to be assessed by both their calibration and discrimination. Calibration addresses the question of how closely the model-based risk estimates align with the observed outcomes. Discrimination focuses on a model’s ability to discriminate individuals who develop the event of interest from those who do not. “Useful” prediction models can then be developed into risk prediction algorithms or rules. Such algorithms can be utilized to classify people into different (e.g., high, intermediate, or low risk; treat pharmacologically, introduce lifestyle intervention, or do not act, respectively) medical decision (or risk) categories based on predetermined thresholds (35). After the marker has been shown to predict incidence of the disease of interest, it must be demonstrated that it adds incremental value to risk prediction models that contain standard factors (generally referred to as the baseline model). Because novel markers must be associated with the incidence of the disease of interest after controlling for risk factors already included in the baseline model, it is necessary to assess their incremental value for risk prediction models in terms other than statistical significance.
The cutpoint-based NRI, cutpoint-free NRI, and IDI each consider separately individuals who do and do not develop events. Therefore, they provide additional information that could not be obtained from discrimination analysis like Harrell’s C index of discrimination. For the NRI, each individual is assigned to a risk category e.g., low (< 6%), medium (6% to < 12%), or high (≥ 12%).
based on the event probability calculated by the baseline model. An enhanced model is developed by adding the biomarker of interest to the reference model, and each individual is reassigned to a risk category. The net proportion of patients with events reassigned to a higher risk category (NRIevents) and of patients without events reassigned to a lower risk category (NRInonevents) is calculated. The NRI is the sum of NRIevents and NRInonevents. It is interpreted as the proportion of patients reclassified to a more appropriate risk category. Among those with the event, if the addition of the biomarker of interest to the model results in more individuals being reclassified to higher risk categories than to lower ones, then the NRIevents is positive. Conversely, among those without events, if more are assigned to lower than higher risk categories, then the NRInonevents is positive (35).
Harrell’s C index of discrimination and Nam-D’Agostino X2 can provide useful information on predictive performance of a predictive model. Examination of the clinical relevance of a new risk biomarker corresponds to the examination of the predictive power of a currently available predictive model augmented by new biomarker(s). NRI and IDI can be extended so that they are applicable to the survival-based regression models. Herein, we have described a Stata package that can be used to calculate these indices. The packages provided herein can encourage novel statistical approaches to be more extensively employed in the study of the predictive capacities of prediction models or the clinical relevance of new biomarkers.
