*! version 2.2.2 06Mar2024

help pmcalplot

Title

pmcalplot -- Calibration plot of prediction model performance

Syntax

pmcalplot varlist [if] [in] [,Pmcalplot_general_options Binary_outcome_options survival Survival_outcome_options continuous Continuous_outcome_options Twoway_options]

where varlist varies dependent on the type of outcome of interest (binary, survival or continuous). See below for requirements and options relevant to each outcome type.
General options available for all outcome types are given below.

pmcalplot_general_options

nostatistics suppresses display of performance statistics

range(# #) a range describing the square size of the plot region

scatteropts(string) twoway options to affect rendition of the groups

nosmoother suppresses display of the smoother

nosmootherci suppresses display of the CI for the smoother

smootheropts(string) twoway options to affect rendition of the smoother

smootherciopts(string) twoway options to affect rendition of the CI for the smoother

keep returns variables for the expected, observed & risk groups

zoom zooms the graph area to fit the risk groups and their CI's

twoway_options all other twoway options

Binary outcomes (Logistic prediction models)

For binary outcome models varlist must consist of two variables; the first variable must represent the predicted probabilities from the model (a variable in the range [0,1]).
The second variable must represent the event indicator i.e. the observed outcome (a binary variable coded 0 for non-events and 1 for events).

binary_outcome_options

bin(int 10) number of bins used to group patients average observed & expected probabilities

cutpoints(numlist) unequal cutpoints used to group patients average observed & expected probabilities

ci displays 95% confidence intervals for groups

ciopts(string) twoway options to affect rendition of the confidence intervals

nospike suppresses display of the spike plot

spikeopts(string) twoway options to affect rendition of the spike plot

Survial outcomes (Cox/RP prediction models)

For survial outcomes varlist must consist of only one variable; the first variable must represent the predicted probabilities from the model (a variable in the range [0,1]).
The observed outcome probabilities are calculated by pmcalplot using Kaplan-Meier estimates.
NOTE: Predicted probabilites provided in varlist must represent the probability of an event, i.e. 1-S(t).

survival_outcome_options

survival option required to tell pmcalplot survival outcomes are being used (uses KM estimates for observed data)

timepoint(int 1) timepoint at which observed & expected probabilities are to be compared

lp(varname numeric) variable name representing the linear predictor values based on the development model, for all patients

bin(int 10) number of bins used to group patients average observed & expected probabilities

cutpoints(numlist) unequal cutpoints used to group patients average observed & expected probabilities

ci displays 95% confidence intervals for groups

ciopts(string) twoway options to affect rendition of the confidence intervals

nospike suppresses display of the spike plot

spikeopts(string) twoway options to affect rendition of the spike plot

ttesmoothergroups(int 5) number of groups used to estimate the smoother for TTE outcomes

Continuous outcomes (Linear prediction models)

For continuous outcomes varlist must consist of two variables; the first variable must represent the predicted outcome values.
The second variable must represent the observed outcome values.

continuous_outcome_options

continuous option required to tell pmcalplot continuous outcomes are being used

p(int 0) number of variables used in development model

nohist suppresses display of the histograms for observed and expected outcomes

obshistopts(string) twoway options to affect rendition of the histogram of observed values

exphistopts(string) twoway options to affect rendition of the histogram of expected values

Description

pmcalplot produces a calibration plot of observed against expected probabilities for assessment of prediction model performance.
Calibration is plotted in groups across the risk spectrum as recommended in the TRIPOD guidelines, and confidence intervals for the groupings can also be displayed (NB: not for continuous outcomes).
Further, a spike plot of the distribution of events and non-events can be displayed on the plot, as well as a smoother allowing assessment of the calibration at the individual patient level.
For continuous outcomes a histogram of observed and expected values can be displayed on the corresponding axes.
Additionally, common prediction model performance statistics can also be displayed on the plot, quantifying the models performance.

pmcalplot is primarily useful for assessment of model performance in an external validation of an existing model.
However it can also be used during model development to check the apparent performance of the model (which should show perfect calibration).

General Options

nostatistics specifies that the model performance statistics are not calculated and displayed on the calibration plot (note this may be computationally faster).

range(# #) specifies a new plot range given two real numbers. The plot range must lie within the bounds of 0 and 1 for binary and survival outcome models. For continuous outcome models the range is dependent on the outcome measure.
Both axes are forced to have the same range so that the plot remains square in line with the TRIPOD guidelines. The default range is 0 to 1 for binary & survival outcome models.

zoom automatically 'zooms' the plot range to fit the risk groups and their CI's. Smoother and spike plots are only plotted within the 'zoomed' ranged. The range() option overides zoom.
Useful for prediction models with a narrow range of predictions.

scatteropts(string) specifies any additional options to describe the rendition of the scatter plot. See scatter.

ciopts(string) specifies any additional options to describe the rendition of the confidence intervals. See twoway_rspike.

nosmoother specifies that the smoother not be calculated and displayed on the calibration plot (note this may be computationally faster for large datasets).

nosmootherci specifies that the 95% confidence interval for the smoother not be calculated and displayed on the calibration plot.

smootheropts(string) specifies any additional options to describe the rendition of the smoother. See twoway_line.

smootherciopts(string) specifies any additional options to describe the rendition of the 95% CI for the smoother. See twoway_rarea.

keep specifies that variables for the expected, observed & risk groups defined by pmcalplot be retained in the dataset after the program terminates.

twoway_options controls all other graph options, e.g. legend layout, axis titles etc. See twoway_options.

Binary outcome options

bin(int 10) specifies the number of equally sized percentiles to divide the predicted risks into. The default is to divide the predicted risks into 10 equally sized groups.

cutpoints(numlist) specifies cutpoints to divide the predicted risks into.
numlist should contain a list of real numbers in the range [0,1] to define risk thresholds for grouping.
The cutpoints() option is a useful alternative to the default when validating an existing prediction model with defined risk thresholds to group patients (allowing unequal group sizes).
It is advisable not to use the cutpoints() option unless the model has prespecified clinically relevant cut-points.

ci specifies that 95% confidence intervals for the observed proportions be added to the plot.

nospike specifies that the spike plot not be displayed on the calibration plot.

spikeopts(string) specifies any additional options to describe the rendition of the spike plot. See twoway_rspike.

Survival outcome options

survival option required to tell pmcalplot survival outcomes are being used.

timepoint(int 1) specifies a single timepoint at which observed & expected probabilities are to be compared. Units of time are defined by the stset command prior to using pmcalplot. Default = 1 unit of time.
NOTE: User must provide predicted probabilities at the specific time point of interest.
This will usually require some estimate of the baseline survival or baseline hazard at the specified time point.
Please consult the TRIPOD guidelines if unsure. See sts_list, stset.

lp(varname numeric) if a variable containing the linear predictor values for patients in the dataset is provided, then some performance statistics can be displayed on the plot.
Such a variable should be available as it is used in the calculation of predicted/expected risks which must be given in varlist for survival outcome models.
The linear predictor variable must be given for the smoother to be calculated and displayed for survival/TTE outcomes.

ttesmoothergroups(int 10) specifies the number of groups used to estimate the smoother via pseudo-values. The default is five groups.
If the dataset is sufficiently large it is recommended to use 10 or more groups.

bin(int 10) specifies the number of equally sized percentiles to divide the predicted risks into. The default is to divide the predicted risks into 10 equally sized groups.

cutpoints(numlist) specifies cutpoints to divide the predicted risks into.
numlist should contain a list of real numbers in the range [0,1] to define risk thresholds for grouping.
The cutpoints() option is a useful alternative to the default when validating an existing prediction model with defined risk thresholds to group patients (allowing unequal group sizes).
It is advisable not to use the cutpoints() option unless the model has prespecified clinically relevant cut-points.

ci specifies that 95% confidence intervals for the observed proportions be added to the plot.

nospike specifies that the spike plot not be displayed on the calibration plot. Spike plot for survival models, plots the events and non-events at the timepoint of interest.

spikeopts(string) specifies any additional options to describe the rendition of the spike plot. See twoway_rspike.

Continuous outcome options

continuous option required to tell pmcalplot continuous outcomes are being used.

p(int 0) if the number of variables included in the development model is provided, then R-squared adjusted can additionally be calculated and displayed on the plot.
NOTE: p is the number of variables in the published model to be validated.

nohist specifies that the histograms not be displayed on the calibration plot.

obshistopts(string) specifies any additional options to describe the rendition of the histogram of observed values (y-axis histogram). See histogram.

exphistopts(string) specifies any additional options to describe the rendition of the histogram of expected values (x-axis histogram). See histogram.

Examples

Binary outcomes (Logistic prediction models)

Load dataset

. webuse lbw, clear

For illustration only, create a binary variable to define a development and validation cohort

. gen val = runiform()<.4

Derive a prediction model for low birth weight in the development cohort

. logistic low age lwt i.race smoke ptl ht ui if val==0

Generate a new variable containing the predicted probabilities from the model for all individuals

. predict p_dev

Use pmcalplot to produce a calibration plot of the apparent performance of the development model in the development cohort

. pmcalplot p_dev low if val==0, ci

Now produce the calibration plot for the models performance in the validation cohort

. pmcalplot p_dev low if val==1, ci

We can also remove elements of the calibration plot such as the spike plot and performance statistics

. pmcalplot p_dev low if val==1, ci nospike nostatistics

Other twoway options can be used to alter the rendition of the calibration plot such as to remove the legend or alter titles

. pmcalplot p_dev low if val==1, ci nospike nostatistics xtitle("Predicted probability", size(medsmall)) ytitle("Observed frequency", size(medsmall)) legend(off)

Now produce the calibration plot for the models performance in the validation cohort using risk thresholds of <5%, 5-15%, 15-50% and >50%. This is useful if the model has been proposed for decision making at these thresholds.
Use the keep option to return the groupings to your dataset

. pmcalplot p_dev low if val==1, cut(.05 .15 .5) ci keep

Survial outcomes (Cox prediction models)

Load dataset

. webuse drugtr, clear

Fit a flexible parametric survival model & predict survival probabilities for individuals at 10 year follow-up

. stset

. stpm2 drug age, df(3) scale(h)

. gen time1 = 10

. predict p_surv, s timevar(time1)

Derive the probability of an event i.e. 1-S(t) (where S(t) is the probability of survival)

. gen p_event = 1-p_surv

Use pmcalplot to produce a calibration plot of the apparent performance of the development model in the development cohort

. pmcalplot p_event, ci surv t(10)

Derive the linear predictor for all individuals

. predict lp, xbnobaseline

Feed the linear predictor to pmcalplot to produce a calibration plot including performance statistics for the apparent performance of the model

. pmcalplot p_event, ci surv t(10) lp(lp)

Continuous outcomes (Linear prediction models)

Load dataset

. sysuse auto, clear

Derive a prediction model for mpg in the development cohort

. regress mpg weight foreign

Generate a new variable containing the predicted values from the model for all individuals

. predict exp

Use pmcalplot to produce a calibration plot of the apparent performance of the development model in the development cohort.
Telling pmcalplot the number of factors in the model also allows R-squared adjusted to be calculated

. pmcalplot exp mpg, cont p(2)

Now set the range & histogram axes to better include all data in pmcalplot

. pmcalplot exp mpg, cont r(10 50) obshi(ylab(10(10)50)) exphi(xlab(10(10)50))

Stored results

pmcalplot stores the following in r() when the option nostatistics is not specified (as relevant for specific outcome types):

Scalars

r(cstat) C-statistic for model discrimination

r(oe_ratio) Ratio of observed and expected events for model calibration

r(citl) Calibration-in-the-large for model calibration

r(cslope) Calibration slope for model calibration

r(r2) R-squared for overall model performance (linear models)

r(r2a) R-squared for overall model performance adjusted for optimism (linear models)

Authors

Joie Ensor, University of Birmingham

j.ensor@bham.ac.uk

Acknowledgements

With thanks to Emma Martin, Kenneth Pihl, Rebecca Whittle & Rupert Major for helpful feedback

References

Ensor J. Calibration plots of prediction model performance in Stata using pmcalplot. In preparation.

Also see

Online: twoway_options, running