%Corrected when received galleys: (Y-Yhat)^2 (added ^2), -.0972->-.0972*age
\documentclass{article}
%\documentstyle[here,citesort,S,feh,psbox]{article}
\usepackage{/home/harrelfe/doc/latex/feh}
\usepackage{/home/harrelfe/doc/latex/texinput/citesort}
\usepackage{/home/harrelfe/doc/latex/texinput/S}
\input{/home/harrelfe/doc/latex/texinput/psbox.tex}
\usepackage{/home/harrelfe/doc/latex/here}
\def\comment#1{\relax}

\newcommand{\ctable}[3]{        	%
        \begin{center}          	%
		\begin{table}[H]			%
		\caption{#2}				%
		\label{#3}					%
        \begin{tabular}{#1}}

\renewcommand{\etable}{				%
		\end{tabular}      			%
		\end{table}					%
        \end{center}            	%
%       \vspace{1.5em}				%
		}

%\input{/tmp/chbars.sty}

\begin{document}
\renewcommand\thefootnote{}	%was \arabic{footnote}
\title{MULTIVARIABLE PROGNOSTIC MODELS:\\ Issues in Developing Models,
Evaluating Assumptions and Adequacy, and Measuring and Reducing Errors\\
\emph{Statistics in Medicine} \textbf{15}:361-387; 1996}
\footnotetext{
From the Division of Biometry, Department of Community and Family Medicine,
and the Division of Cardiology, Department of Medicine, Duke University
Medical Center.

This work was supported by: Research Grants HL-17670, HL-29436,
HL-36587, HL-45702 and HL-09315 from the National Heart, Lung and
Blood Institute, Bethesda, Maryland; Research Grants HS-03834,
HS-05635, HS-06503, HS-06830, and HS-07137 from the Agency for Health
Care Policy and Research, Rockville, Maryland; and grants from the
Robert Wood Johnson Foundation, Princeton, NJ.

Address for correspondence and reprint requests:
\begin{tabbing}
~~~~~ \= Frank E Harrell Jr, PhD \\
\> Box 3363 \\
\> Duke University Medical Center \\
\> Durham, North Carolina 27710 USA\\
\> E-mail: \Co{feh@biostat.mc.duke.edu} \\
\end{tabbing}
}
\renewcommand\thefootnote{\arabic{footnote}}

\author{Frank E. Harrell, Jr., Ph.D.
\and Kerry L. Lee, Ph.D. \and Daniel B. Mark, M.D., M.P.H.}
\date{\today}
\maketitle
\noindent{\bf ABSTRACT}
\vspace{1ex} \\

Multivariable regression models are powerful tools that are used
frequently in studies of clinical outcomes. These models can
use a mixture of categorical and continuous variables and can handle
partially observed (censored) responses. However, uncritical
application of modeling techniques can result in models that poorly fit
the dataset at hand or, even more likely, inaccurately predict outcomes
on new subjects.

One must know how to measure qualities of a model's fit in order to
avoid poorly fitted or over--fitted models. Measurement of predictive
accuracy can be difficult for survival time data in the presence of
censoring. We discuss an easily interpretable index of predictive
discrimination as well as methods for assessing calibration of
predicted survival probabilities. Both types of predictive accuracy
should be unbiasedly validated using bootstrapping or
cross--validation, before using predictions in a new data series.

We discuss some of the hazards of poorly fitted and overfitted
regression models and present one modeling strategy that avoids many
of the problems discussed. The methods described are applicable to all
regression models, but are particularly needed for binary, ordinal, and
time--to--event outcomes.  Methods are illustrated with a survival
analysis in prostate cancer using Cox regression.
\begin{tabbing}
\underline{Key words}: \=
prognosis, prediction, multivariable regression models, survival
analysis \\
\> predictive accuracy, validation, calibration, bootstrap, modeling strategy
\end{tabbing}

\section{Introduction} 
Accurate estimation of patient prognosis is important for many
reasons.  First, prognostic estimates can be used to inform the
patient about likely outcomes of her disease. Second, the physician
can use estimates of prognosis as a guide for ordering additional
tests and selecting appropriate therapies. Third, prognostic
assessments are useful in the evaluation of technologies; prognostic
estimates derived both with and without using the results of a given
test can be compared to measure the incremental prognostic information
provided by that test over what is provided by prior information
\cite{har82}. Fourth, a researcher may want to estimate the effect of
a single factor (e.g., treatment given) on prognosis in an
observational study in which many uncontrolled confounding factors are
also measured. Here the simultaneous effects of the uncontrolled
variables must be controlled (held constant mathematically if using a
regression model) so that the effect of the factor of interest can be
more purely estimated. An analysis of how variables (especially
continuous ones) affect the patient outcomes of interest is necessary
to ascertain how to control their effects.  Fifth, prognostic
estimation is useful in designing randomized clinical trials. Both the
decision concerning which patients to randomize and the design of the
randomization process (e.g., stratified randomization using prognostic
factors) are aided by the availability of accurate prognostic
estimates before randomization \cite{spi86}.  Lastly, accurate
prognostic models can be used to test for differential therapeutic
benefit or to estimate the clinical benefit for an individual patient
in a clinical trial, taking into account the fact that low--risk
patients must have less absolute benefit (lower change in survival
probability) \cite{kna93cli}.

To accomplish these objectives, analysts must create prognostic models
that accurately reflect the patterns existing in the underlying data
and that are valid when applied to comparable data in other settings
or institutions.  Models may be inaccurate due to violation of
assumptions, omission of important predictors, high frequency of
missing data and/or improper imputation methods, and especially with
small datasets, overfitting. The purpose of this paper is to review
methods for examining lack of fit and detection of overfitting of
models and to suggest guidelines for maximizing model accuracy.
Section 2 covers initial steps such as imputation of missing data,
pre--specification of interactions, and choosing the outcome model.
Section 3 has an overview of the need for data reduction. In Section
4, we discuss the process of checking whether a hypothesized model
fits the data. In Section 5, measures of predictive accuracy are
covered. These are not directly related to lack of fit but rather to
the ability of the model to discriminate and be well calibrated when
applied prospectively.  Section 6 covers model validation and
demonstrates advantages of resampling techniques.  Section 7 provides
one modeling strategy that takes into account
ideas from earlier sections and lists some miscellaneous concerns.
Most of the methods presented here can be used with any regression
model.  
Section 8 briefly describes some statistical software useful in
carrying out the strategy summarized in Section 7.  Section 9
has a detailed case study using a Cox regression
model for time until death in a clinical trial studying prostate
cancer.

\section{Preliminary Steps}
Before analyses begin, the researcher must specify the relationships
of interest and define and assemble the response variable and the
potential predictors.  At this point a frequent problem is the extent
of missing data.  Some methods of dealing with missing data are given
in \cite{don82rel,rob87sas,buc60met,tim70est}.  Deletion of cases with
missing predictors causes bias and increased variance.  Even though
caution should be taken when imputing missing values, it is usually
better to estimate selected data values than to delete an entire
subject's record.  Simple methods of imputation include the use of the
median, mean, or mode for missing values.  This method is biased and
inefficient when predictors are correlated with one another
\cite{don82rel}.  Deriving customized regression models for predicting
each predictor from all other predictors is a better method.  Kuhfeld
\cite{prinqual} has implemented a general imputation method that
allows predictors to be nonlinearly (and even non--monotonically)
related to one another.  This method has been modified by Harrell and
implemented in the S--Plus {\tt transcan} function (Section
\ref{software}), which yields stable imputations even when the
fraction of missing values is quite large.  In some cases, surrogate
predictors, not intended to enter the model directly, are assembled to
assist in imputing missing predictors in the model.

It is important that maximum information be extracted from predictors
and response.  Because of this and because of problems with data
reliability, when one has a choice of describing a concept with a
categorical variable or a continuous one, the continuous one is preferred.
Subject matter knowledge should guide the selection of candidate
predictors.  Early deletion of those with little chance of being
predictive or of being measured reliably will result in models with
less overfitting and greater generalizability.  

Plausible interactions
should be carefully chosen because of problems of multiple parameters
(see \cite{sch88non} for additional thoughts on interactions).
Certain types of interactions that have frequently been found to be
important in predicting clinical outcomes and thus may be
pre--specified are:
\begin{enumerate}
\item Interactions between treatment and the severity of disease being treated.
Patients with little disease have little opportunity to receive benefit.
\item Interactions involving age and risk factors.  Old subjects are generally
less affected by risk factors.  They have been robust enough to survive to
their current age with risk factors present.
\item Interactions involving age and type of disease.	Some diseases are
incurable and have the same prognosis regardless of age.  Others are
treatable or have less effect on younger patients.
\item Interactions between a measurement and the state of a subject during
a measurement.  For example, left ventricular function measured at
rest may have less predictive value and thus have a smaller slope
vs. outcome than function measured during stress.
\item Interactions between calendar time and treatment.  Some treatments
evolve or their effectiveness improves with staff training.
\item Interactions between quality and quantity of a symptom.
\end{enumerate}
Careful fitting of a statistical model is essential so that
interactions, if present, represent biologic phenomena rather than
general lack--of--fit of the model.

A tentative choice of the statistical model is sometimes based on
previous distributional examinations, but it is frequently based on
maximizing how available information is used.  Binary and ordinal
logistic models \cite{cox58reg,wal67,hou88log,col91} are frequently
used for discrete completely assessed outcomes, and the Cox
proportional hazards model \cite{cox72,col94} and parametric survival
models \cite{law82} are frequently used for censored time--to--event
data.  It is quite common to change the model after initial modeling
of predictors, because only then can adjusted distributional
properties of $Y$ and joint properties of $X$ and $Y$ be assessed
(Section \ref{distributional}).

\section{Data Reduction}
Multivariable statistical models when developed carefully are
excellent tools for making prognostic predictions. However, when the
assumptions of a model are grossly violated or when a model is used
unwisely for a given patient sample, the performance of the model may
be poor. For example, when the analyst has fitted not only real trends
that further data would support, but in addition has fitted
idiosyncrasies in the particular dataset by analyzing too many
variables, the model may predict inaccurately for a new group of
patients. Only with appropriate model validation can an apparently
accurate model be shown to be inaccurate.

In developing a set of predictions based on 100 patients, no analyst
would divide the patients into 50 subgroups and quote the average
outcome for each subgroup. Yet many articles have appeared in the
clinical literature where 20--50 variables were analyzed on 100
patients. Researchers apparently do not realize that when many
predictor variables are analyzed, variable screening based on
statistical significance and stepwise variable selection involve
multiple comparisons problems that lead to unreliable models.  These
methods are therefore not viable for data reduction. [See
\cite{der92bac} for a condemnation of stepwise variable selection.]

The situation is actually worse than merely considering the number of
predictors.  If the analyst used associations with $Y$ to entertain
nonlinearities in the predictors or interaction terms, these
constructed variables need to be counted (see Table \ref{predictors}
for an example).  We speak of the total predictor degrees of freedom
(d.f.\ ), $p$, as the total number of parameters (columns of the
design matrix) examined during the course of analysis, excluding
intercept term(s).  If graphical or other informal analyses are used
to guide the analysis, it is difficult to define $p$ --- one needs to
estimate the effective number of parameters considered according to
the flexibility of fits that were considered \cite{gra91}.  The
quantity $p$ is the effective number of parameters allowed for
consideration, i.e., the number of regression coefficients estimated
formally or informally without algebraic restrictions.  

To enhance the accuracy of a model, the number of variables used must
be reduced or the model must be simplified unless the sample is large.
Unless a formal penalized estimation technique is used \cite{ver94pen},
multiple comparisons problems that arise from ``peeking'' at
the outcome variable must be eliminated; data reduction methods must
be used that do not utilize the outcome variable. Harrell et
al. \cite{har84} discussed some available data reduction methods and
two regression modeling strategies based on these methods that yield
reliable models. They suggest as a rough rule of thumb that 
in order to have predictive discrimination that validates on a new
sample, no more than $m/10$ predictor d.f.\ $p$
should be examined to fit a multiple
regression model, where $m$ is the number of uncensored event times
(e.g., deaths) in the training sample 
(the sample used in fitting the model).
For binary outcomes $m$ is the number of patients in
the less frequent outcome category. If $p > m/10$,
a data reduction technique such as principal
components, variable clustering, or deriving clinical summary indexes
\cite{mar94ass,har84,jol86,jac91use}
should be used until the number of summary variables to use as
candidates in the regression analysis is less than $m/10$.

Smith et al. \cite{smi92pro} found in one series of simulations that
the expected error
\footnote{Absolute difference between predicted and actual 5--year survival
probabilities in a simulation study with known survival functions}
in Cox model predicted 5--year survival
probabilities was below 0.05 when $p < m/20$ for ``average'' subjects
and below 0.10 when $p < m/20$ for ``sick'' subjects.  For ``average'' subjects, $m/10$ was adequate for preventing
expected errors $>0.1$.

Better and more general than any of these rules is the reduction of
d.f.\ using a shrinkage method (Section \ref{shrinkage}).


\section{Verifying Model Assumptions: Checking Lack of Fit}

\subsection{Linearity Assumption}
In their simplest forms, all usual regression models assume that for a
certain scale of $Y$, each predictor variable $X$ is linearly related to
$Y$.  In the logistic regression model for binary responses,
the initial assumption is that an $X$
is linearly related to the log odds of response ($\log[P/(1-P)]$, where
$P$ is the probability of response) for patients subgrouped by values
of $X$.  In the Cox proportional hazards survival model,
one initially assumes that at each time $t$, $\log[-\log(S(t))]$ and
equivalently $\log \lambda(t)$ are linearly related to $X$, where
$S(t)$ is the probability of surviving until time $t$ and $\lambda(t)$
is the hazard function or instantaneous event rate at time $t$.	It is
easy to envision cases where strong violations in the linearity
assumption (say a U--shaped age relationship) will result in erroneous
predictions.

A direct way to check the linearity assumption, and to determine how to
transform a specific $X$ if necessary, involves expanding $X$ into
multiple terms that can flexibly fit any smooth relationship. The extra
terms can be statistically tested to assess the adequacy of a linear
relationship, and the terms in toto can estimate the true
transformation of $X$ that would result in a linear relationship with
$Y$. A common choice of expansion is to add $X^2$ and perhaps higher
powers of $X$ to the model. A more flexible approach is the use of
piecewise linear regression or piecewise cubic polynomials (spline
functions). See \cite{dur89,har88,sle90} for methods of fitting such
functions.

As an alternative, smoothed residual plots  can be used to
determine the functional form for each predictor.  For binary logistic
models, smoothed partial residual plots
\cite{lan84gra,col91,hos89} are useful, and for the Cox model,
smoothed martingale residuals plots detect regression shape departures
\cite{the90}.	Partial residuals in logistic models are particularly
computationally efficient, as the analyst can fit a simple model
that is linear in all predictors and then use the residuals to
obtain estimates of the {\em true} functional forms.  However, the
plot for each predictor does assume that the other predictors operate
linearly and that all predictors are additive (see below).
The usual martingale residual plot for the Cox model provides an
estimate of the {\em departure} from linearity for the predictor.

\subsection{Additivity Assumption}
\label{additivity}
A further assumption of most regression models is additivity of
effects of the predictors (lack of interaction). Interactions can
be tested and described by adding cross--product terms. It must
be borne in mind that interactions can take the form of a change
in shape (e.g., linear age relationship for males, quadratic for
females), so the cross--products needed in the model are not
always simple ones.

The number of possible cross--product terms is usually so large
(especially when variables have nonlinear or multiple dummy variable
components)
that the predictors to check for additivity must usually be
specified before examining the data. Otherwise, type I errors
and overfitting will be significant problems.  A  compromise solution
is to do pooled interaction tests.  For example, in a model with
predictors age, sex, and dose, one may test all second--order interactions
involving age, all interactions involving sex, and all involving dose.
A combined test of all
two--way interactions is also useful.  If a pooled test is not significant,
it may be unwise to pursue significant component interactions.

\subsection{Distributional Assumption}
\label{distributional}
The previous sections dealt with the proper specification of the $X$--
structure of the model.
Once the analyst has determined which predictors are to be used
and how they should be represented in the model, most models have
distributional assumptions that also need verification. The Cox
model does not assume anything about the survival
function $S(t)$ across $t$ for an individual, but it does assume
how survival curves for different subjects are related. Specifically,
it assumes that $\log[-\log(S(t))]$ for different subjects are
equidistant over time, or equivalently that hazard functions for
any two subjects are proportional over time. This proportional hazards
assumption can be checked using smoothed plots of a special type
of residual from the model called the Schoenfeld residual
\cite{sch82,pet90inv}.
It can also be checked using hazard ratio plots, plots of modeled
vs. stratified estimates
\footnote{That is, a Cox model is fitted with the variable in question
appearing as a covariate for which regression coefficient(s) are
estimated, then a second model is fitted where that variable is used
as a stratification factor that modifies the underlying survival
function (but which does not have regression coefficients).}
, and several other methods \cite{har87}.
Unlike the Cox model, fully parametric models (e.g., Weibull or
log--normal survival models) have a distributional
assumption even when there are no covariables. If the form of $S(t)$
does not fit the data for these models, estimates of $S(t)$ will be
inaccurate.


\section{Quantifying Predictive Accuracy}

There are at least three uses of measures of predictive accuracy:
\begin{enumerate}
\item 
To quantify the utility of a predictor or model to be used for prediction
or for screening to identify subjects at increased risk of a disease
or clinical outcome.
\footnote{Often one wishes to designate a model as
``minimally acceptable'' on the basis of some statistic, but in many
cases it is only possible to judge a model's accuracy relative to
another model. For example, a model for the probability of death
after open heart surgery may yield predicted probabilities that
range from 0.001 to 0.1, so the model will not have a high correlation
(say 0.13) between predicted probability and observed outcome, but it may still
be useful.  If that model does not fully adjust for patient risk
factors, it may be inadequate for adjusting for case mix
when comparing mortalities among several hospitals.  A more sensitive
model with a correlation of, say, 0.135 may adjust away apparent
differences in mortality among hospitals.}
\item To check a given model for overfitting (fitting noise resulting
in unstable regression coefficients) or lack of fit (improper model
specification, omitted predictors, or underfitting).  More will
be said about this later.
\item To rank competing methods or competing models.
\end{enumerate}

The measures discussed below may be applied to the assessment of
a predictive model using the same sample on which the model was
developed.  However, this assessment is seldom of interest, as
only the most serious lack of fit will make a model appear not to
fit on the sample for which it was tailor--made. Of much greater
value is the assessment of accuracy on a separate sample or a bias--corrected
estimate of accuracy on the training sample. This
assessment can detect gross lack of fit as well as overfitting,
whereas the {\em apparent} accuracy from the original model
development sample does not allow one to quantify overfitting.
Section \ref{validation} discusses how the indexes described
below may be estimated fairly using a validation technique.

\subsection{General Notions}
In the simplest case, when the response
being predicted is a continuous variable that is measured completely
(as distinct from {\em censored} measurements caused by termination of
follow--up before all subjects have had the outcome of interest), one
commonly used measure of predictive accuracy is the {\em expected squared error}
of the estimate. This quantity is defined as the expected squared
difference between predicted and observed values, that is, the average
squared difference between predicted and observed values if the
experiment were repeated infinitely often and new estimates were made
at each replication. The expected squared error can also be expressed
as the square of the {\em bias} of the estimate plus the {\em variance}
of the estimate.  Here  bias refers to the expected value of the
estimate minus the quantity being estimated, such as the
mean blood pressure.  The expected squared error is estimated in
practice by the usual mean squared error.

There are two other terms for describing the components of
predictive accuracy: {\em calibration} and {\em discrimination}.
Calibration refers to the extent of bias. For example, if the average
predicted mortality for a group of similar patients is 0.3 and
the actual proportion dying is 0.3, the predictions are well
calibrated.  Discrimination measures a predictor's ability to
separate patients with different responses. A weather forecaster
who predicts a 0.15 chance of rain every day of the year may be
well calibrated in a certain locality if the average number of
days with rain is 55 per year, but the forecasts are
uninformative.  A discriminating forecaster would be one who
assigns a wide distribution of predictions and whose predicted
risks for days where rain actually occurred are larger than for dry
days.  If a predictive model has poor discrimination, no
adjustment or calibration can correct the model. However, if
discrimination is good, the predictor can be calibrated without
sacrificing the discrimination [see section \ref{validation} for
a method for calibrating predictions without needing more data].
Here, calibrating predictions means modifying them, without
changing their rank--order, such that the predictions are
perfectly calibrated. van Houwelingen and le Cessie
\cite{hou90pre} present extensive information on predictive
accuracy and model validation.

\subsection{Continuous Uncensored Outcomes}
Discrimination is related to the expected squared error and to the
correlation between predicted and observed responses. In the case of
ordinary multiple linear regression, discrimination can be measured by the
squared multiple correlation coefficient $R^{2}$, which is defined by
\beq
R^{2} = 1-(n-p){\rm MSE}/(n-1)S_{Y}^{2},
\eeq
where $n$ is the number of patients, $p$ is the number of
parameters estimated, MSE is the mean squared error of prediction
($\sum_{i=1}^{n} (Y_{i}-\hat{Y}_{i})^{2}/(n-p), \hat{Y}$= predicted
$Y$), and $S_{Y}^{2}$ is the sample variance of the dependent
variable. When $R^{2} = 1$, the model is perfectly able to
separate all patient responses based on the predictor variables,
and MSE$=0$.

For a continuous uncensored response $Y$, calibration can be assessed by a
scatter plot of $\hat{Y}$ (predicted $Y$) vs. $Y$, optionally
using a nonparametric smoother to make trends more evident.

\subsection{Discrete or Censored Outcomes}
When the outcome variable is dichotomous and predictions are
stated as probabilities that an event will occur, calibration and
discrimination are more informative than expected squared error
alone in measuring accuracy. 

One way to assess calibration of probability predictions is to form
subgroups of patients and check for bias by comparing predicted and observed
responses \cite[pp. 140-145]{hos89}.  For example, one may group by deciles
of predicted probabilities and plot the mean response (proportion with the
outcome) versus the mean prediction in the decile group. However, the
groupings can be quite arbitrary. Another approach is to use a smoother such
as the ``super smoother'' \cite{fri84} or a scatterplot smoother
\cite{cle79} to obtain a nonparametric estimate of the relationship
between $\hat{Y}$ and $Y$. Such smoothers work well even when $Y$ is
binary. The resulting smoothed function is a nonparametric calibration
or reliability curve. Smoothers operate on the raw data ($\hat{Y}, Y$) and do not require grouping $\hat{Y}$, but
they do require one to choose a smoothing parameter or bandwidth.

As an example, consider a 7--variable binary
logistic regression model to predict the probability that a certain disease is
present. The model was developed on a simulated 200--subject dataset of whom 93
had a final diagnosis that is positive.  While fixing the intercept and 7
regression coefficients estimated from the training sample, predictive
probabilities of disease were computed for each of 200 subjects in a
separate sample, of whom 104 had the disease. The nonparametric
calibration curve was estimated using a local least squares scatterplot smoother \cite{cle79} with the
S--Plus function {\tt lowess} \cite{Splus}, using the ``no iteration''
option. The smoothed calibration graph is shown in Figure \ref{cal.graph}.
Also shown is the proportion of patients with disease,
grouped by intervals of predicted probability each containing 50
patients.
\fig{cal.graph}{Smooth nonparametric calibration curve (dashed line),
subgroup estimates (dots), and ideal relationship (dotted line).
The distribution of predicted probabilities is shown above the
$x$--axis.
``Actual Probability'' is an unbiased estimate of the true probability of response given the level of the predicted probability.}
Note the typical regression to the mean effect caused by overfitting:
predicted probabilities in the range of 0.3--0.5 are too low. Actual
probabilities are closer to the mean (104/200 = 0.52).

When $Y$ is binary and $\hat{Y}$ is the predicted probability that
$Y=1$ vs. $Y=0$, the Brier score \cite{bri50ver} or average
$(Y-\hat{Y})^{2}$ is a commonly used mean squared error--type measure
of predictive accuracy.

For survival models, one may choose one or more times
($t_{1},t_{2},\ldots,t_{k}$), and plot the
predicted probability of surviving until each $t_{j}$ vs. the actual
fraction of patients surviving past $t_{j}$.  The problem here is
that we cannot define $Y_{i}=1$ if patient $i$ survives past time $t_{j}$
and then plot the mean $Y$ (by deciles of $\hat{Y}$ or using a smoother)
against the mean $\hat{Y}$, since subjects not followed until time $t_{j}$ are
censored, i.e., their final outcome status is unknown.  One solution
is to divide the sample into intervals of $\hat{Y}$ so that there are
50 subjects in each interval of predicted survival, and then plot the
mean $\hat{Y}$ within each interval vs. the Kaplan--Meier \cite{kap58}
survival estimate at time $t_{j}$.



\subsection{Shrinkage}
\label{shrinkage}
{\em Shrinkage} is the flattening of the plot of (predicted, observed)
away from the $45^{\circ}$ line, caused by overfitting.  It is a
concept related to regression to the mean.  One can estimate the
amount of shrinkage present (using external validation) or the amount
likely to be present (using bootstrapping, cross--validation, or
simple heuristics).  A shrinkage coefficient can be used to quantify
overfitting or one can go a step further and use the coefficient to 
re--calibrate the model.
Shrinkage can be defined as a multiplier 
$\gamma$ of $X\hat{\beta}$ (excluding
intercept(s)) needed to make $\gamma X\hat{\beta}$ perfectly calibrated
for future data. The heuristic shrinkage estimator of van Houwelingen
and le Cessie \cite{hou90pre} (see also \cite{cop83reg}) is 
\beq
\hat{\gamma} = \frac{{\rm model}\ \chi^{2} - p}{{\rm model}\ \chi^{2}},
\label{heuristic.shrink}
\eeq
where $p$ is the number of regression parameters (here excluding any
intercept(s) but including all nonlinear and interaction effects) and
the model $\chi^2$ is the total likelihood ratio $\chi^2$ statistic
(computed using the full set of $p$ parameters)
for testing whether any predictors are associated with $Y$
\footnote{When stepwise fitting is done, the definition of $p$ is
confusing.  Many analysts act as if the final model chosen with
stepwise variable selection was pre--specified, whether interpreting
$R^2$, confidence limits, or $P$--values.  For estimating the
likely shrinkage, it has been shown that $p$ is much closer to the number of
candidate d.f.\ than to the number of parameters fitted
in a ``final'' model \cite{cop83reg}.  On a similar note, \cite{gra91}
showed how to adjust a linear test of association for having done a
test of quadratic effect, concluding that testing the single d.f.\ 
statistic for association as if it had 2 d.f.\ is nearly optimal.}
. For linear regression, van
Houwelingen and le Cessie's heuristic shrinkage estimate reduces to
the ratio of the adjusted $R^2$ to the ordinary $R^2$ 
(derivable from \cite[Eq. 70]{hou90pre}).

As an example, suppose that an analyst has considered 10
predictor variables,
6 of which were allowed to enter the model non--linearly
(with 2 nonlinear terms for each), and tested 8 interaction terms, for
a total of 30 degrees of freedom.  The model $\chi^2$ is 100 for the
full model fit with $p=30$ d.f. 
The expected shrinkage is 0.70, indicating that about
0.3 of the model fit is ``noise''. The ``final model'' obtained from
forward variable selection contains only 3
significant coefficients and has $\chi^2 = 81$, but overfitting is
quantified using the 30 candidate d.f.  In this example, the number
of variables, transformations, and interactions tried was too many for
the sample size, and the resulting model is expected to be unstable.
As a rough estimate, 0.3 of what was learned from developing the model
was really non--replicable noise.

For mild overfitting in the case where the model is needed only to
rank likely outcomes and not predict absolute risks, shrinking the
regression coefficients will not help since it won't increase real
discrimination. If the model is badly overfitted, the model may actually have
negative (worse than random) discrimination on new data, and it will
have poor calibration. The following heuristic strategy can then be
used to determine whether data reduction is likely to result in a
model that has any discrimination and how much reduction is required
to yield reliable non--shrunken predictions.

First, fit a full model with all candidate variables, nonlinear
terms, and hypothesized interactions. Let $p$ denote the number
of parameters in this model, aside from any intercept(s). Let LR
denote the likelihood ratio $\chi^{2}$ for this full model.
The estimated shrinkage is $({\rm LR} - p)/{\rm LR}$. If this
falls below 0.85, for example, we may be concerned.  Let $q$
denote the regression degrees of freedom for a reduced model. In
a ``best case'', the variables removed to arrive at the reduced
model would have no association with $Y$. The expected value of
the $\chi^{2}$ statistic for testing those variables would then
be $p-q$. The shrinkage for the reduced model is then on average
$[{\rm LR} - (p-q) - q]/[{\rm LR} - (p-q)]$.  Solving for $q$
gives $q \leq ({\rm LR}-p)/9$. Therefore, reduction of
dimensionality down to $q$ degrees of freedom would be expected
to achieve $<10\%$ shrinkage. With these assumptions, there is no
hope that a reduced model would have acceptable calibration
unless LR $> p+9$. If the information explained by the omitted
variables is less than one would expect by chance (e.g., their
total $\chi^{2}$ is extremely small), a reduced model could still
be beneficial, as long as the conservative bound $({\rm LR} -
q)/{\rm LR} \geq 0.9$ or $q \leq {\rm LR}/10$ were achieved. This
conservative bound assumes that no $\chi^{2}$ is lost by the
reduction, i.e., that the final model $\chi^{2} \approx$ LR. This
is unlikely in practice, since the data reduction {\em must} be
only $X$--driven.

As an example, suppose that a binary logistic model is being developed
from a sample containing 45 events on 150 subjects. A 10:1
events:d.f.\  rule
suggests we can analyze 4.5 degrees of freedom. The analyst wishes to
analyze age, sex, and 10 other variables. It is not known whether
interaction between age and sex exists, and whether age is linear.
A restricted cubic spline is fitted with 4 knots (requiring two nonlinear
terms), and a linear interaction
is allowed between age and sex. These two variables then need $3+1+1 = 5$
degrees of freedom. The other 10 variables are assumed to be linear
and to not interact with themselves or age and sex. There is a total of
15 d.f.  The full model with 15 d.f.\  has LR=50. Expected shrinkage
from this model is $(50-15)/50 = 0.7$.	Since LR$>15+9=24$, some
reduction {\em might} yield a better validating model. Reduction to
$q = (50-15)/9 \approx 4$ d.f.\  would be necessary, assuming the reduced
LR is about $50-(15-4) = 39$. In this case the 10:1 rule yields
about the same value for $q$. The analyst may be forced to assume
that age is linear, modeling 3 d.f.\  for age and sex. The other 10 variables
would have to be reduced to a single variable using principal components
or another scaling technique. This single variable may not be
interpretable, but using a single score is better than deleting all 10
variables from consideration.  If the goal of the analysis is to make
a series of hypothesis tests (adjusting $P$--values for multiple
comparisons) instead of to predict future responses, the full model
would have to be used.

Bootstrapping \cite{hou90pre} and cross--validation \cite{cop87cro}
may also be used to estimate shrinkage factors.

As mentioned above, shrinkage estimates are useful in their own right
for quantifying overfitting, and they are also useful for ``tilting''
the predictions so that the (predicted, observed) plot does follow the
$45^{\circ}$ line, by multiplying all of the regression coefficients
by $\hat{\gamma}$.  However, for the latter use it is better to follow
a more rigorous approach such as penalized maximum
likelihood estimation \cite{ver94pen}, which allows the analyst to
shrink different parts (e.g., nonlinear terms or interactions) of the
equation more than other parts \cite{gra92fle}.

\subsection{General Discrimination Index}

Discrimination can be defined more uniquely than calibration.  It can be
quantified with a measure of correlation without requiring the formation of
subgroups or requiring smoothing.

When dealing with binary dependent variables or continuous dependent
variables that may be censored when some patients have not suffered the
event of interest, the usual mean squared error--type measures do not apply.
A $c$ (for {\em concordance}) index \cite{har82} is a
widely applicable measure of predictive discrimination---one that applies
to ordinary continuous outcomes, dichotomous diagnostic outcomes, ordinal
outcomes, and censored time--until event response variables. This index of
predictive discrimination is related to a rank correlation between
predicted and observed outcomes. It is a modification of the
Kendall--Goodman--Kruskal--Somers type rank correlation index \cite{goo79} and was
motivated by a modification of Kendall's $\tau$ by Brown, Hollander, and
Korwar \cite{bro74non} and Schemper \cite{sch84ana}.

The $c$ index is defined as the proportion of all usable patient pairs in
which the predictions and outcomes are concordant. The $c$ index measures
predictive information derived from a set of predictor variables in a
model. In predicting the time until death, $c$ is calculated by considering
all possible pairs of patients, at least one of whom has died. If the
predicted survival time is larger for the patient who lived longer, the
predictions for that pair are said to be concordant with the outcomes. If
one patient died and the other is known to have survived at least to the
survival time of the first, the second patient is assumed to outlive the
first. When predicted survivals are identical for a patient pair, $\frac{1}{2}$
rather than 1 is added to the count of concordant pairs in the numerator of
$c$. In this case, one is still added to the denominator of $c$ (such patient
pairs are still considered usable). A patient pair is unusable if both
patients died at the same time, or if one died and the other is still alive
but has not been followed long enough to determine whether she will outlive
the one who died.

Instead of using the predicted survival time to calculate $c$, the predicted
probability of surviving until any fixed time point can be used
equivalently, as long as the two estimates are one--to--one functions of each
other. This holds for example if the proportional hazards assumption is
satisfied.

For predicting binary outcomes such as the presence of disease, $c$ reduces
to the proportion of all pairs of patients, one with and one without the
disease, in which the patient having the disease had the higher predicted
probability of disease. As before, pairs of patients having the same
predicted probability get $\frac{1}{2}$ added to the numerator. The
denominator is the number of patients with disease multiplied by the number
without disease. In this binary outcome case, $c$ is essentially the
Wilcoxon--Mann--Whitney statistic for comparing predictions in the two
outcome groups, and it is identical to the area under a receiver operating
characteristic (ROC) curve \cite{bam75,han82}. Liu and Dyer \cite{liu79ran}
advocate the use of rank association measures such as $c$ in quantifying the
impact of risk factors in epidemiologic studies.

The $c$ index estimates the probability of concordance between predicted and
observed responses. A value of 0.5 indicates no predictive discrimination
and a value of 1.0 indicates perfect separation of patients with different
outcomes. For those who prefer instead a rank correlation coefficient
ranging from $-1$ to $+1$ with 0 indicating no correlation, Somers' $D$
rank correlation index is derived by calculating $2(c-0.5)$. Either $c$ or
the rank correlation index can be used to quantify the predictive
discrimination of any quantitative predictive method, whether the response
is continuous, ordinal, or binary.

Even though rank indexes such as $c$ are widely applicable and
easily interpretable, they are not sensitive for detecting small
differences in discrimination ability between two models. This is
due to the fact that a rank method considers the (prediction,
outcome) pairs $(0.01, 0), (0.9, 1)$ as no more concordant than
the pairs $(0.05, 0), (0.8, 1)$. A more sensitive
likelihood--ratio $\chi^{2}$ -- based statistic that reduces to
$R^2$ in the linear regression case may be substituted
\cite{nag91not,lee86,har85com}. Korn and Simon \cite{kor90mea}
have a very nice discussion of various indexes of accuracy for
survival models.

\section{Model Validation Methods}
\label{validation}

As mentioned before, examination of the {\em apparent} accuracy of a
multivariable model using the training dataset is not very useful.
The most stringent test of a model (and of the entire data collection
system) is an external validation --- the application of the
``frozen'' model to a new population.  It is often the case that the
failure of a model to validate externally could have been predicted
from an honest (unbiased) ``internal'' validation.  In other words, it is
likely that many clinical models which failed to validate would have
been found to fail on another series of subjects from the original
source, because overfitting is such a common problem.  The principal
methods for obtaining nearly unbiased internal assessments of accuracy
are {\em data--splitting} \cite{pic90}, {\em cross--validation}
\cite{efr83est}, and {\em bootstrapping}
\cite{efr83est,lin89ass,efr83lei,efr86,boot}.  In data--splitting, a
random portion, for example $\frac{2}{3}$, of the sample is used for
all model development (data transformations, stepwise variable
selection, testing interactions, estimating regression coefficients,
etc.).  That model is ``frozen'' and applied to the remaining sample
for computing calibration statistics, $c$, etc.  The size of the
validation sample must be such that the relationship between predicted
and observed outcomes can be estimated with good accuracy, and the
remaining data are used as the training 
(model development)
sample.  Data--splitting is
simple, because all the modeling steps, which may include subjective
judgements, are only done once.  Data--splitting also has an advantage
when it is feasible to make the single split with respect to
geographical location or time, resulting in a more stringent
validation that demonstrates generalizability.  However,
in addition to severe difficulties listed below, data splitting does
not validate the final model, if one desires to recombine the training
and test data
to derive a model for others to use.

Cross--validation is repeated data--splitting.  To obtain accurate
estimates using cross--validation, more than 200 models may need to be
developed and tested \cite{efr83est}, with results averaged over the 200
repetitions.  For example, in a sample of size
$n=1000$, the modeling process (all components of it!) could be done
400 times, leaving out a random 50 subjects each time and developing
the model on the 950 remaining subjects.  The benefits of
cross--validation over data--splitting are clear: the size of the
training samples can be much larger, so less data are discarded from
the estimation process; secondly, cross--validation reduces variability
by not relying on a single sample split.

Efron has shown that cross--validation is relatively inefficient due
to high variation of accuracy estimates when the entire validation
process is repeated \cite{efr83est}. Data--splitting is far worse; the
indexes of accuracy will vary greatly with different
splits. Bootstrapping is an alternative method of internal validation
that involves taking a large number of samples with replacement from
the original sample. Bootstrapping provides nearly unbiased estimates
of predictive accuracy that are of relatively low variance, and fewer
model fits are required than cross--validation. Bootstrapping has an
additional advantage that the entire dataset is used for model
development.  As others have shown, data are too precious to waste
\cite{roe91pre,bre92lit}.  

Suppose that we wish to estimate the expected value (for new patient
samples similar to the derivation sample) of the Somers' $D$ rank
correlation coefficient between predicted and observed survival
time. The following steps can be used (see
\cite{lin89ass,bre92lit,boot} for the basic method when applied to
binary outcomes):
\begin{enumerate}
\item Develop the model using all $n$ subjects and whatever stepwise testing is deemed necessary. Let $D_{\rm app}$ denote the {\em apparent} $D$ from this model, i.e., the rank correlation computed on the same sample used to derive the fit.
\item Generate a sample of size $n$ with replacement from the original sample (for both predictors and the response).
\label{step1}
\item Fit the full or possibly stepwise model, using the same stopping rule as was used to derive $D_{\rm app}$.
\item Compute the apparent $D$ for this model on the bootstrap sample with replacement. Call it $D_{\rm boot}$.
\item ``Freeze'' this reduced model, and evaluate its performance on the original dataset. Let $D_{\rm orig}$ denote the $D$.
\item The optimism in the fit from the bootstrap sample is $D_{\rm boot} - D_{\rm orig}$.
\label{step2}
\item Repeat steps \ref{step1} through \ref{step2} 100--200 times.
\item Average the optimism estimates to arrive at $O$.
\item The bootstrap corrected performance of the original
stepwise model is $D_{\rm app} - O$. This difference is a nearly unbiased
estimate of the {\em expected value} of the external predictive
discrimination of the {process} which generated $D_{\rm app}$.
In other words, $D_{\rm app} - O$ is an {\em honest} estimate
of {\em internal} validity, penalizing for overfitting.
\end{enumerate}
As an example, suppose we want to validate a stepwise Cox model
developed from, say, a sample of size $n=300$ with 30 events. The
candidate regressors are age, age$^2$, sex, mean arterial blood
pressure (MBP), and a nonlinear interaction between age and sex with
the terms age $\times$ sex and age$^2 \times$ sex. MBP is assumed to
be linear and additive.  Denote these variables by the numbers
1--6. The model $\chi^2$ is 45 with 6 d.f.\ , so the approximate
expected shrinkage is $\frac{45-6}{45} = 0.87$, or 0.13 overfitting,
so some caution needs to be exercised in using the estimated model
coefficients and hence in using extreme predicted survival
probabilities without calibration (shrinkage).  The $D$ for the full model is
0.42. A step--down variable selection using Akaike's information
criterion (AIC) \cite{atk80,hou90pre} as a stopping rule ($\chi^2$ for
set of variables tested $> 2 \ \times$ d.f.\ ) resulted in a model with
the variables age, age$^2$, sex, age $\times$ sex. The reduced model
had $D=0.39$, a typical loss due to deleting marginally important but
statistically insignificant variables. Two--hundred bootstrap
repetitions are done, repeating the variable selection for each sample
using the same stopping rule. We want to detect whether the $D=0.39$
is likely to validate in a new series of subjects from the same
population. The first 5 samples might yield the following results:

\ctable{cccccc}{Example validation of predictive discrimination}{discrimval} 
\hline\hline
Re--sample & $D_{\rm boot}$ & Variables Retained & $D_{\rm boot}$ &
$D_{orig}$ & Optimism \\ & Full Model & & Reduced Model & & \\ \hline
1 & 0.45 & 1,2,3,5,6 & 0.44 & 0.37 & 0.07 \\ 2 & 0.46 & 1,2 & 0.34 &
0.30 & 0.04 \\ 3 & 0.42 & 1,2,3,4 & 0.37 & 0.34 & 0.03 \\ 4 & 0.43 &
1,2,3,5 & 0.42 & 0.39 & 0.03 \\ 5 & 0.41 & 1,3,4 & 0.39 & 0.37 & 0.02
\\ \hline \etable The average optimism is 0.038, so the bootstrap
estimate of the expected validation of $D_{\rm app}$ is $0.39 - 0.038
= 0.352$.  The analyst may or may not be worried about the 0.038
overfitting, but the best estimate of predictive discrimination
is $D = 0.352$ --- this is a better estimate of the likely
``external'' validation accuracy than is 0.39 if all other aspects of
the study design remain constant.  The $D=0.352$ is the honest estimate of
predictive accuracy that should be quoted when the researchers
document the accuracy of the reduced model that was developed
on the entire dataset using a stepwise variable selection algorithm.

If is usually informative to repeat the bootstrap validation with
and without stepwise variable selection. Usually, the amount of
predictive information lost by deleting marginal variables is not
offset by the decreased optimism of the stepwise model.  One way to
demonstrate this point is to observe how often ``insignificant''
clinical predictors have clinically sensible signs on their
regression coefficients.  Stepwise variable selection, which requires
binary decisions about the inclusion of variables (unlike shrinkage),
causes information to be lost \cite{spi86}.

The same strategy can be used to estimate the over--optimism in
an $R^{2}$ measure \cite{nag91not} from the original model fit.
For estimating the prediction error at time $t$ in a survival
model, similar steps could also be used. Instead of validating a correlation
$D$, we substitute for example the statistic $D$ =
difference between mean predicted 2--year survival probability
and Kaplan--Meier 2--year survival estimate. The survival
estimates are made by, say, deciles of predicted 2--year survival from
the original model fit using the following steps, for example.

\begin{enumerate}
\item Develop the model using all subjects
\item Compute cut points on predicted survival at 2 years so that there
	are $m$ patients within each interval ($m=50$ or $100$ typically)
\item For each interval of predicted probability, compute the mean
	predicted 2--year survival and the Kaplan--Meier 2--year
	survival estimate for the group
\item Save the apparent errors --- the differences between mean predicted
	and Kaplan--Meier survival
\item Generate a sample with replacement from the original sample
	\label{step.a}
\item Fit the full model
\item Do variable selection and fit the reduced model
\item Predict 2--year survival probability for each subject in the bootstrap
	sample
\item Stratify predictions into intervals using the previously chosen cut points
\item Compute Kaplan--Meier survival at 2 years for each interval
\item Compute the difference between the mean predicted survival within
	each interval and the Kaplan--Meier estimate for the interval
\item Predict 2--year survival probability for each subject in the
	original sample using the model developed on the
	sample with replacement
\item For the same cut points used before, compute the difference in the
	mean predicted 2--year survival and the corresponding
	Kaplan--Meier estimates for each group in the original sample
\item Compute the differences in the differences between the bootstrap
	sample and the original sample
	\label{step.b}
\item Repeat steps \ref{step.a} to \ref{step.b} 100--200 times
\item Average the ``double differences'' computed in step \ref{step.b}
	over the 100--200 bootstrap samples.
	These are the estimates of over--optimism in the apparent error
	estimates.
\item Add these over--optimism estimates to the apparent errors in the
	original sample to obtain
	bias--corrected estimates of predicted vs. observed, i.e., to
	obtain a bias-- or overfitting--corrected calibration curve.
\end{enumerate}

\section{Summary of Modeling Strategy}
\label{strategy}
\begin{enumerate}
\item	Assemble accurate, pertinent data and as large a sample as possible.
	For survival
	time data, follow--up must be sufficient to capture enough
	events as well as the clinically meaningful phases if dealing
	with a chronic disease.
\item	Formulate focused clinical hypotheses that lead to specification of
	relevant candidate predictors, the form of expected relationships,
	and possible interactions.
\item	Discard observations having missing $Y$ 
	after characterizing
	whether they are missing at random. \footnote{For survival time data,
	no observations should be missing on $Y$. They should only
	have curtailed follow--up.}
See \cite{cra95com} for a study of imputation of $Y$ when it is not
missing at random.
\item	If there are any missing $X$s, analyze factors associated with
	missing--ness. If the fraction of observations that would be
	excluded due to missing values is very small, or one of the variables
	that is sometimes missing is of overriding importance, exclude
	observations with missing values\footnote{Alternatively, impute
	missing values for the predictor but perform secondary analyses
	later to estimate the strength of association between $X$ and $Y$
	after deleting observations with that predictor imputed, as imputation will
	attenuate the relationship.}. Otherwise impute missing $X$s using
	individual predictive models that take into account the reasons
	for missing, to the extent possible.
\item	If the number of terms fitted \underline{or} tested in the
	modeling process (counting nonlinear and cross--product terms)
	is too large in comparison with the number of outcomes in
	the sample, use data reduction (ignoring $Y$) 
\cite{mar94ass,har84,jol86,jac91use}
	until the
	number of remaining free variables needing regression coefficients
	is tolerable. Assessment of likely shrinkage (overfitting) can be
	useful in deciding how much data reduction is adequate.
	Alternatively, build shrinkage into the initial model fitting
	\cite{ver94pen}.
\item	Use the entire sample in the model development as data are too
	precious to waste. If steps listed below are too difficult to
	repeat for each bootstrap or cross--validation sample, hold
	out test data from all model development steps which follow.
\item	Check linearity assumptions and make transformations in $X$s
	as needed. \label{step.beg}
\item	Check additivity assumptions and add clinically motivated
	interaction terms.
\item	Check to see if there are overly--influential observations
	\cite{the90}.	Such observations may indicate overfitting,
	the need for truncating the range of highly skewed variables
	or making other pre--fitting transformations,
	or the presence of data errors.
\item	Check distributional assumptions and choose a different model
	if needed (in the case of Cox models, stratification or
	time--dependent covariables can be used if proportional
	hazards is violated).
\item	Do limited backwards step--down variable selection
	\cite{man70why}.  Note that since stepwise techniques do not
	really address overfitting and they can result in a loss of
	information, full model fits (i.e., leaving all hypothesized
	variables in the model regardless of $P$--values) are frequently
	more discriminating than fits after screening predictors for
	significance \cite{cop83reg,spi86}.  They also provide confidence
	intervals with the proper coverage, unlike models that are reduced
	using a stepwise procedure
	\cite{alt89,hur90,bre92lit}, from which confidence
	intervals are falsely narrow.  A compromise would be to test a
	{\em pre--specified} subset of predictors, deleting them if their
	total $\chi^{2} < 2 \ \times\ $ d.f.  If the $\chi^{2}$ is that small, the
	subset would likely not improve model accuracy.
\label{step.end}
\item	This is the ``final'' model.
\item	Validate this model for calibration and discrimination ability,
	preferably using bootstrapping. Steps \ref{step.beg} to
	\ref{step.end} must be repeated for each bootstrap sample,
	at least approximately. For example, if age was transformed
	when building the final model, and the transformation was
	suggested by the data using a fit involving age and age$^{2}$,
	each bootstrap repetition should include both age variables
	with a possible step--down from the quadratic to the linear
	model based on automatic significance testing at each step.
\item	If doing stepwise variable selection, present a summary table
	depicting the variability of the list of ``important factors''
	selected over the bootstrap samples or cross--validations.
	This is an excellent tool for understanding why data--driven
	variable selection is inherently ambiguous.
\item	Estimate the likely shrinkage of predictions from the model,
	either using Eq. \ref{heuristic.shrink} or by bootstrapping
	an overall slope correction for the predictions \cite{hou90pre}.
	Consider shrinking the predictions to make them calibrate better,
	unless shrinkage was built--in.
	That way, a predicted 0.4 mortality is more likely to validate
	in a new patient series, instead of finding that the actual mortality is
	only 0.2 because of regression to the mean mortality of 0.1.
\end{enumerate}

\section{Software}
\label{software}
Modern statistical software such as S-Plus \cite{Splus} on UNIX
workstations makes it quite feasible to perform the extensive
calculations required to do the recommended model building steps.  
The first author has written a package of UNIX S--Plus functions
called {\tt Design} \cite{Design} that allow the analyst to perform all
analyses mentioned here including tests of linearity, pooled
interaction tests, model validation, and graphical methods for
interpreting models.  Here are some examples:
\need{.6in}
\begin{Example}
    # First find optimum transformations relating each predictor to each
    # other, and use multiple regression in these transformations to
    # impute missing values.  Use shrinkage to avoid over-imputing
    trans \Gets transcan({\Twiddle} age + cholesterol + sys.bp + weight, imputed=T, shrink=T)
    cholesterol \Gets impute(trans, cholesterol)  # impute missings
    sys.bp \Gets impute(trans, sys.bp)
    # Fit a Cox P.H. model allowing some interactions with age and
    # nonlinearity in cholesterol and sys.bp using restricted cubic splines
    # x=T, y=T means store data in fit for future bootstrapping
    fit \Gets cph( Surv(fu.time, death) {\Twiddle} age * (rcs(cholesterol) + rcs(sys.bp)) +
                         weight, x=T, y=T, surv=T, time.inc=5)
    anova(fit)                      # automatic pooled Wald tests
    fastbw(fit)                     # fast backward step-down
    # Next validate model, penalizing for backward stepdown variable selection
    validate(fit, B=100, bw=T)      # bootstrap validation of accuracy indexes
    calibrate(fit, B=100, bw=T, u=5)# bias-corrected 5-yr survival calibration
    plot(summary(fit))              # plot hazard ratios with confidence limits
    nomogram(fit)                   # draw nomogram displaying how model works
    latex(fit)                      # typeset model equation
\end{Example}
The {\tt Design} library includes a function {\tt rcorr.cens} for
computing the general $c$--index, and the function {\tt val.prob}
which produced Figure \ref{cal.graph} and also prints a variety of accuracy
measures.  For binary and ordinal logistic models and for ordinary
linear models, {\tt Design} has a general penalized maximum likelihood
estimation facility.
{\tt Design} is available in the {\tt statlib} repository ({\tt
Internet} address {\tt lib.stat.cmu.edu}).  {\tt transcan} and
\Co{impute} are separate functions in {\tt statlib} which work on
UNIX as well as DOS Windows S--Plus.  Some other software systems
which have some intermediate--level capabilities include Stata
(Computer Resources Center Inc., College Station TX), SPIDA (NHMRC
Clinical Trials Centre, Eastwood, NSW Australia), and SAS (SAS
Institute Inc., Cary NC).

\section{Case Study}
Consider the 506--patient prostate cancer dataset from Byar and Green
\cite{bya80cho} which has also been analyzed in
\cite{kay86tre,sau92boo}.  The data are listed in \cite[Table
46]{data} and are available by Internet at \Co{utstat.toronto.edu} in
the directory \Co{/pub/data-collect}. These data were from a randomized
trial comparing four treatments for stage 3 and 4 prostate cancer,
with almost equal numbers of patients on placebo and each of three
doses of estrogen.  Four patients had missing values on all of the
following variables: \Co{wt, pf, hx, sbp, dbp, ekg, hg, bm}; two of
these patients were also missing \Co{sz}
(see Table \ref{predictors} for abbreviations).
These patients will be excluded from consideration.

There are 354 deaths among the 502 patients.  If we only wanted to
test for a drug effect on survival time, a simple rank--based analysis
would suffice.  To be able to test for differential treatment effect
or to estimate prognosis or expected absolute treatment benefit for
individual patients, however, we need a multivariable survival model 
\cite{kna93cli}.
First we consider fitting a full additive model which
does not assume linearity of effect for any predictor. Categorical
predictors will be expanded using dummy variables.  For \Co{pf} we
could lump the last two categories since the last category has only 2
patients. Likewise, we could combine the last two levels of \Co{ekg}.
Continuous predictors will be expanded by fitting 4--knot restricted
cubic spline functions, which contain 2 nonlinear terms and thus have
a total of 3 d.f.   The following table defines the candidate
predictors and lists their d.f.  The variable \Co{stage} is not listed
as it can be predicted with high accuracy from \Co{sz,sg,ap,bm}
(\Co{stage} could have been used as a predictor for imputing
missing values on \Co{sz, sg}).

\ctable{llcl}{Candidate predictors and d.f.}{predictors} \hline\hline
\multicolumn{1}{c}{Predictor} & \multicolumn{1}{c}{Name} & Number of &
Original Levels\\
 & & Parameters & \\ \hline
Dose of estrogen					&	\Co{rx}		& 3 	& 
placebo, 0.2, 1.0, 5.0 mg estrogen							  \\
Age in years						&	\Co{age}	& 3 	& \\
Weight index: wt(kg)--ht(cm)+200	&	\Co{wt}		& 3 	& \\
Performance rating					&	\Co{pf}		& 2 	& 
	normal, in bed $<50\%$ of time, 						  \\
& & & in bed $>50\%$, in bed always							  \\
History of cardiovascular disease	&	\Co{hx}		& 1 	&present/absent\\
Systolic blood pressure/10			&	\Co{sbp}	& 3 	& \\
Diastolic blood pressure/10			&	\Co{dbp}	& 3 	& \\
Electrocardiogram code				&	\Co{ekg}	& 5 	&
	normal, benign, rhythm disturb., block, 				  \\
& & & strain, old myocardial infarct, new MI 				  \\
Serum hemoglobin (g/100ml)			&	\Co{hg}		& 3 	& \\
Tumor size (cm$^2$)					&	\Co{sz}		& 3 	& \\
Stage/histologic grade combination	&	\Co{sg}		& 3		& \\
Serum prostatic acid phosphatase	&	\Co{ap}		& 3		& \\
Bone metastasis						&	\Co{bm}		& 1		&present/absent\\
\hline
\etable

There are a total of 36 candidate d.f.\  which should not be
artificially reduced by ``univariable screening'' or graphical
assessments of association with death.  This is about $\frac{1}{10}$
as many predictor d.f.\  as there are deaths, so there is some hope that
a fitted model may validate.  Let us also examine this issue by
estimating the amount of shrinkage using Equation
\ref{heuristic.shrink}.  We use a Cox proportional hazards model for
time until death.  The UNIX S--Plus \Co{Design} library fits the
full model using restricted cubic spline expansions and makes use of
Therneau's \Co{survival4} package in \Co{statlib} \cite{survival} to perform the calculations.
But first we invoke the \Co{transcan} function and \Co{impute} functions
(from \Co{statlib} for any versions of S--Plus) to develop customized
nonlinear imputation equations for all predictors and to apply these
equations to impute missing values.
\begin{Example}
    #Define function for easy determination of whether a value is in a list
    '%in%' \Gets function(a,b) match(a, b, nomatch=0) > 0

    levels(ekg)[levels(ekg) %in% c('old MI','recent MI')] \Gets 'MI'
    # combines last 2 levels and uses a new name, MI

    pf.coded \Gets as.integer(pf)    # save original pf, re-code to 1-4
    levels(pf)  \Gets c(levels(pf)[1:3], levels(pf)[3]) # combine last 2 levels of original

    w \Gets transcan({\Twiddle} sz + sg + ap + sbp + dbp + age + wt + hg + 
                    ekg + pf + bm + hx, imputed=T, impcat='tree')
    sz  \Gets impute(w, sz)          # uses imputation rule w
    sg  \Gets impute(w, sg)
    age \Gets impute(w, age)
    wt  \Gets impute(w, wt)
    ekg \Gets impute(w, ekg)

    dd \Gets datadist(rx, age, wt, pf, pf.coded, heart, map, hg, sz, sg, ap, bm)
    options(datadist='dd')    # datadist stores characteristics of raw data

    units(dtime) \Gets 'Month'
    S \Gets Surv(dtime, status!='alive')

    f \Gets cph(S {\Twiddle} rx + rcs(age,4) + rcs(wt,4) + pf + hx + 
             rcs(sbp,4) + rcs(dbp,4) + ekg + rcs(hg,4) + 
             rcs(sg,4) + rcs(sz,4) + rcs(ap,4) + bm)
\end{Example}
The likelihood ratio $\chi^2$ statistic is 140 with 36 d.f.  This test
is highly significant so some modeling is warranted.  The AIC value
(on the $\chi^2$ scale) is $140 - 2 \times 36 = 68$.  The rough
shrinkage estimate is 0.743 (104/140) so we estimate that 26\% of the
model fitting will be noise, especially with regard to calibration
accuracy.  The approach of \cite{spi86} is to fit this full model and
to shrink predicted values.  We will instead try to do data reduction
(blinded to individual $\chi^2$ statistics from the above model fit)
to see if a reliable model can be obtained without shrinkage.  A good
approach at this point might be to perform a variable clustering analysis
which for our purposes we will do informally.  The data reduction
strategy is listed in Table \ref{reduce.strategy}.  For \Co{ap}, more
exploration is desired to be able to model the shape of effect with
such a highly skewed distribution. Since we expect the tumor variables
to be strong prognostic factors we will retain them as separate
variables.  No assumption will be made for the dose--response shape
for estrogen, as there was reason to expect a non--monotonic effect
due to competing risks for cardiovascular death.

\ctable{llr}{Data reduction strategy (blinded to $Y$)}{reduce.strategy} 
\hline\hline
\multicolumn{1}{c}{Variables} & 
\multicolumn{1}{c}{Reductions} & 
\multicolumn{1}{c}{d.f.\  Saved} \\ \hline
\Co{wt}		&	Assume variable not important enough for 4 knots & 1 \\
			&	Use 3 knots										 &   \\
\Co{pf}		&	Assume linearity								 & 1 \\
\Co{hx,ekg}	&	Make new 0,1,2 variable and assume linearity:	 & 5 \\
			&	2=\Co{hx} and \Co{ekg} not normal or benign,	 &	\\
			&   1=either, 0=none			 					 & \\
\Co{sbp,dbp}&	Combine into mean arterial bp and use 3 knots:	 & 4 \\
			&	\Co{map}=$\frac{2}{3}$ \Co{dbp} $+ \frac{1}{3}$ \Co{sbp} & \\
\Co{sg}		&	Use 3 knots										 & 1 \\
\Co{sz}		&	Use 3 knots										 & 1 \\
\Co{ap}		&	Look at shape of effect of \Co{ap} in detail,	 & -2 \\
			& 	and take log before expanding in spline to achieve & \\
			& 	numerical stability: add 2 knots				 & \\ \hline
\etable

\need 2in
\begin{Example}
    heart \Gets hx + !( ekg %in% c('normal','benign') )
    label(heart) \Gets 'Heart Disease Code'
    map   \Gets (2*dbp + sbp)/3
    label(map) \Gets 'Mean Arterial Pressure/10'

    f \Gets cph(S {\Twiddle} rx + rcs(age,4) + rcs(wt,3) + pf.coded + 
             heart + rcs(map,3) + rcs(hg,4) + 
             rcs(sg,3) + rcs(sz,3) + rcs(log(ap),6) + bm,
             x=T, y=T, surv=T, time.inc=5*12)
    # x, y for predict, validate, calibrate;  surv, time.inc for calibrate
\end{Example}

The total savings is thus 11 d.f.  The likelihood ratio $\chi^2$ is
126 with 25 d.f.\ , with a slightly improved AIC of $76$.  The rough
shrinkage estimate is slightly better at $0.80$, but still worrisome.
A further data reduction might be achieved by using the \Co{transcan}
transformations determined from self--consistency of predictors, but
we will stop here and use this model.

\comment{Of the four tumor variables \Co{sg} and \Co{ap} can be predicted from
the other variables (after appropriate transformations) with $R^2$ of
$0.56$ and $0.57$, respectively, and the first 3 principal components
(PCs) of the four transformed variables explain 0.97 of the variance.
The first PC explains 0.65 of the total variance.  Let us see how much
predictive accuracy is lost by summarizing the four variables into a
single score using the first PC.  The PC is computed on transformed
variables (\Co{transcan} used a method similar to that used in
\cite{prinqual} which maximizes the proportion of variance explained
by the first PC).}

Now assess this model in more detail by examining coefficients and
summarizing multiple parameters within predictors using Wald statistics.

\begin{Example}
    f        # writing an object name in S causes it to be printed

Cox Proportional Hazards Model

cph(formula = S {\Twiddle} rx + rcs(age, 4) + rcs(wt, 3) + pf.coded + heart + rcs(map, 3) +
        rcs(hg, 4) + rcs(sz, 3) + rcs(sg, 3) + rcs(log(ap), 6) + bm,
        x=T, y=T, surv=T, time.inc=5*12)

 Obs Events Model L.R. d.f. P Score Score P    R2 
 502    354        126   25 0   135       0 0.221


                        coef se(coef)       z        p 
rx=0.2 mg estrogen  3.74e-03 1.50e-01  0.0250 9.80e-01
rx=1.0 mg estrogen -4.21e-01 1.66e-01 -2.5427 1.10e-02
rx=5.0 mg estrogen -9.73e-02 1.58e-01 -0.6176 5.37e-01
               age -1.17e-02 2.35e-02 -0.4995 6.17e-01
              age'  2.00e-02 3.86e-02  0.5190 6.04e-01
             age''  2.71e-01 4.95e-01  0.5482 5.84e-01
                wt -2.46e-02 9.39e-03 -2.6175 8.86e-03
               wt'  1.84e-02 1.12e-02  1.6379 1.01e-01
          pf.coded  2.25e-01 1.21e-01  1.8625 6.25e-02
             heart  4.18e-01 8.08e-02  5.1723 2.31e-07
               map  3.24e-02 8.49e-02  0.3817 7.03e-01
              map' -4.57e-02 9.41e-02 -0.4857 6.27e-01
                hg -1.56e-01 7.68e-02 -2.0343 4.19e-02
               hg'  7.42e-02 2.10e-01  0.3530 7.24e-01
              hg''  5.08e-01 1.27e+00  0.4014 6.88e-01
                sz  1.00e-02 1.44e-02  0.6955 4.87e-01
               sz'  8.79e-03 2.37e-02  0.3715 7.10e-01
                sg  7.19e-02 7.86e-02  0.9138 3.61e-01
               sg' -7.04e-03 9.83e-02 -0.0716 9.43e-01
                ap -7.96e-01 3.11e-01 -2.5584 1.05e-02
               ap'  4.89e+01 2.18e+01  2.2482 2.46e-02
              ap'' -3.64e+02 1.59e+02 -2.2909 2.20e-02
             ap'''  4.04e+02 1.75e+02  2.3057 2.11e-02
            ap'''' -9.69e+01 4.16e+01 -2.3311 1.97e-02
                bm  3.25e-02 1.81e-01  0.1790 8.58e-01

    # The terms with ', '', etc. after the name are cubic spline nonlinear terms
    # The dose effect is apparently nonlinear.

    anova(f)    # output was actually typesetted automatically using latex(anova(f))
                # latex requires the print.display package from statlib
\end{Example}

% latex.default(x = dstats, title = title, rowlabel = "", caption = head, cdec =  
%
\begin{table}[H]
 \begin{center}
 \caption{Wald Statistics for {\tt S }\label{anova.f}} 
 \begin{tabular}{lrrr} \hline
\multicolumn{1}{l}{}&
\multicolumn{1}{c}{$\chi^2$}&
\multicolumn{1}{c}{$d.f.$}&
\multicolumn{1}{c}{$P$}
\\ \hline
rx&$  8.38$&$ 3$&~0.0387  \\
age&$ 12.85$&$ 3$&~0.0050  \\
~~{\it Nonlinear}&$  8.18$&$ 2$&~0.0168  \\
wt&$  8.87$&$ 2$&~0.0118  \\
~~{\it Nonlinear}&$  2.68$&$ 1$&~0.1014  \\
pf.coded&$  3.47$&$ 1$&~0.0625  \\
heart&$ 26.75$&$ 1$&$<0.0001$\\
map&$  0.25$&$ 2$&~0.8803  \\
~~{\it Nonlinear}&$  0.24$&$ 1$&~0.6272  \\
hg&$ 11.85$&$ 3$&~0.0079  \\
~~{\it Nonlinear}&$  6.92$&$ 2$&~0.0314  \\
sz&$ 10.60$&$ 2$&~0.0050  \\
~~{\it Nonlinear}&$  0.14$&$ 1$&~0.7102  \\
sg&$  3.14$&$ 2$&~0.2082  \\
~~{\it Nonlinear}&$  0.01$&$ 1$&~0.9429  \\
ap&$ 13.17$&$ 5$&~0.0218  \\
~~{\it Nonlinear}&$ 12.93$&$ 4$&~0.0116  \\
bm&$  0.03$&$ 1$&~0.8579  \\
{\bf TOTAL NONLINEAR }&$ 30.28$&$12$&~0.0025  \\
{\bf TOTAL }&$128.08$&$25$&$<0.0001$\\
\hline
\end{tabular}
\end{center}

\end{table}

There are 12 parameters associated with nonlinear effects, and the
overall test of linearity indicates the strong presence of
nonlinearity for at least one of the variables
\Co{age,wt,map,hg,sz,sg,ap}.  There is a difference in survival time between 
at least two of the doses of estrogen.  

Now that we have a tentative model, let us examine the model's
distributional assumptions.  As mentioned in Section
\ref{distributional}, the Schoenfeld partial residuals are an
effective tool for checking the proportional hazards assumption in the
Cox model.  Grambsch and Therneau \cite{gra94pro} have modified these
residuals so that smoothed plots of them estimate the effect of
predictors on the log instantaneous hazard rate as a function of
follow--up time.  Their scaled residuals estimate $\beta(t)$, the
regression coefficient as a function of time.  A messy detail is how
to handle multiple regression coefficients per predictor.  Here we do
an approximate analysis in which each predictor is scored by adding up
all the terms in the model to transform that predictor to be optimally
related to the log hazard (at least if the {\em shape} of the effect
does not change with time).  In doing this we are temporarily ignoring
the fact that the individual regression coefficients were estimated
from the data.  For dose of estrogen, for example, we code the effect
as 0 (placebo), 0.0037 (0.2 mg), -0.421 (1.0 mg), and -0.0973 (5.0
mg), and \Co{age} is transformed as -.0117 \Co{age} + 0.02 \Co{age'}
+0.271 \Co{age''}, which in most simple form is

\need 0.5in
\begin{eqnarray*}
&  -1.17\!\times\!10^{-2}{\tt age}+ 3.48\!\times\!10^{-5}({\tt age}-56)_{+}^{3}+ 4.71\!\times\!10^{-4}({\tt age}-71)_{+}^{3}  \\
&   -1.01\!\times\!10^{-3}({\tt age}-75)_{+}^{3}+ 5.09\!\times\!10^{-4}({\tt age}-80)_{+}^{3}  \\
\end{eqnarray*}
where $(x)_{+}$ means to ignore that term if $x \leq 0$, and the knots
for
\Co{age} are 56, 71, 75, and 80 years.

In S--Plus the \Co{predict} function easily summarizes multiple terms
and produces a matrix (here, \Co{z}) containing the total effects for
each predictor.  Matrix factors can easily be included in model formulas.
\begin{Example}
    z \Gets predict(f, type='terms')     # required x=T above to store design matrix
    f.short \Gets cph(S {\Twiddle} z, x=T, y=T)  # store x, y so can get residuals
\end{Example}
The fit \Co{f.short} based on the matrix \Co{z} of single d.f.\  predictors has the same LR $\chi^2$ of 126 as the fit \Co{f}, but with a falsely low
11 d.f.  All regression coefficients are unity.

Now get scaled Schoenfeld residuals separately for each predictor and
test the proportional hazards assumption for each using the
``correlation with time'' test.  Also plot smoothed trends in the
residuals.  The \Co{plot} method for \Co{cox.zph} objects uses
restricted cubic splines to smooth the relationship.
\begin{Example}
    phtest \Gets cox.zph(f.short, transform='identity')
    phtest

              rho   chisq      p 
      rx  0.12965  6.5451 0.0105
     age -0.08911  2.8518 0.0913
      wt -0.00878  0.0269 0.8697
pf.coded -0.06238  1.4278 0.2321
   heart  0.01017  0.0451 0.8319
     map  0.03928  0.4998 0.4796
      hg -0.06678  1.7368 0.1876
      sz -0.05262  0.9834 0.3214
      sg -0.04276  0.6474 0.4210
      ap  0.01237  0.0558 0.8133
      bm  0.04891  0.9241 0.3364
  GLOBAL       NA 15.3776 0.1659
\end{Example}
Only the drug effect significantly changes over time ($P=0.01$ for
testing the correlation \Co{rho} between the scaled Schoenfeld
residual and time), but when a global test of PH is done penalizing
for 11 d.f., the $P$ value is 0.17.  A graphical examination of the
trends doesn't find anything interesting for the last 10 variables.  A
residual plot is drawn for \Co{rx} alone and is shown in Figure \ref{rx.ph}.
\begin{Example}
    plot(phtest, var='rx')
\end{Example}

\fig{rx.ph}{Raw and spline--smoothed scaled Schoenfeld residuals for dose of 
estrogen, nonlinearly coded from the Cox model fit, with $\pm$ 2
standard errors \protect{\cite{survival}}.}

We will ignore the possible increase in effect of estrogen over time.
If this non--PH is real, a more accurate model might be obtained by
stratifying on \Co{rx} or by using a time $\times$ \Co{rx} interaction
as a time--dependent covariable.

Note that the model has several insignificant predictors.  These will
not be deleted, as that would not improve predictive accuracy and it
would make confidence intervals for $\hat{\beta}$ or for
predicted survival probabilities with the correct coverage
probabilities hard to obtain \cite{alt89}.  At this point
it would be reasonable to test pre--specified interactions.  Here
we will test all interactions with dose.  Since the multiple terms for
many of the predictors (and for \Co{rx}) make for a great number of d.f.\ 
for testing interaction (and a loss of power), we will do approximate
tests on the data--driven codings of predictors.  $P$--values for these
tests are likely to be somewhat anti--conservative.

\begin{Example}
    z.dose \Gets z[,"rx"]  # same as saying z[,1] - get first column
    z.other \Gets z[,-1]   # all but the first column of z
    f.ia \Gets cph(S {\Twiddle} z.dose * z.other)
    anova(f.ia)

                                      Factor     Chi-Square d.f.     P 
 z.dose  (Factor+Higher Order Factors)            18.9      11   0.062
  All Interactions                                12.2      10   0.273
 z.other  (Factor+Higher Order Factors)          134.3      20   0.000
  All Interactions                                12.2      10   0.273
 z.dose * z.other  (Factor+Higher Order Factors)  12.2      10   0.273
 TOTAL                                           137.3      21   0.000
\end{Example}
Here ``Factor+Higher Order Factors'' means the combined main effect
and interaction effect.
The global test of additivity has $P=0.27$, so we will ignore the
interactions (and also forget to penalize for having looked for them below!).

The following UNIX S--Plus statements plot how each predictor is
related to the log hazard of death, along with 0.95 confidence bands.
Note that due to a peculiarity of the Cox model the standard error of
the predicted $X\hat{\beta}$ is zero at the reference values (medians
here, for continuous predictors).

\begin{Example}
    par(mfrow=c(3,4))    # 4 x 3 matrix of graphs
    r <- c(-1,1)         # use common y-axis range for all
    plot(f, rx=NA,       ylim=r)    # NA -> use default range for predictor
    plot(f, age=NA,      ylim=r)
    scat1d(age)                     # scat1d from statlib, for any S-Plus
    plot(f, wt=NA,       ylim=r)    # scat1d shows data density
    . . .
\end{Example}

\fig{cox.shapes}{Shape of each predictor on log hazard of death.  $Y$--axis
shows $X\hat{\beta}$, but the predictors not plotted are set to
reference values.  ``Rug plots'' on the top of each graph show the
data density of the predictor.  Note the highly non--monotonic
relationship with \protect{\Co{ap}}, and the increased slope after age 70 which
has been found in outcome models for various diseases.}


We first validate this model for Somers' $D_{xy}$ rank correlation
between predicted log hazard and observed survival time, and for slope
shrinkage.  The bootstrap is used (with 200 re--samples) to penalize
for possible overfitting, as discussed in Section \ref{validation}.
\need 1.2in
\begin{Example}
    validate(f, B=200, dxy=T, pr=T)

      index.orig  training     test optimism index.corrected   n 
  Dxy  -0.337377 -0.364644 -0.30976 -0.05488        -0.28250 200
   R2   0.221444  0.261369  0.18445  0.07691         0.14453 200
Slope   1.000000  1.000000  0.78464  0.21536         0.78464 200
\end{Example}

Here ``training'' refers to accuracy when evaluated on the bootstrap
sample used to fit the model, and ``test'' refers to the accuracy when
this model is applied without modification to the original sample.
The apparent $D_{xy}$ is -0.34, but a better estimate of how well the
model will discriminate prognoses in the future is $D_{xy}=-0.28$.
The bootstrap estimate of slope shrinkage is 0.78, surprisingly close
to the simple heuristic estimate.  The shrinkage coefficient could
easily be used to shrink predictions to yield better calibration.

Finally, we validate the model (without using the shrinkage coefficient)
for calibration accuracy in predicting the probability of surviving 5
years.  As detailed in Section \ref{validation}, the bootstrap is used
to estimate the optimism in how well predicted 5--year survival from the
final Cox model tracks Kaplan--Meier 5--year estimates, stratifying by
grouping patients in subsets with about 70 patients per interval of
predicted 5--year survival.

\begin{Example}
    plot(calibrate(f, B=200, u=5*12, m=70))
\end{Example}
%
%
%calibrate.cph(fit = f, u = 5 * 12, m = 70, B = 200)
%
%     index.orig training    test mean.optimism mean.corrected   n 
%[1,]   -0.00749  0.01663  0.0577     -0.037003        0.02951 185
%[2,]    0.01729 -0.00567  0.0109     -0.029964        0.04725 199
%[3,]   -0.04261 -0.00367  0.0219     -0.020202       -0.02240 200
%[4,]   -0.00965 -0.04743 -0.0320      0.000298       -0.00995 200
%[5,]   -0.02419 -0.01553 -0.0264      0.008957       -0.03314 200
%[6,]    0.02834  0.04143 -0.0169      0.032414       -0.00407 200
%[7,]    0.02784  0.01456 -0.0254      0.047881       -0.02004 200
%     mean.predicted     KM KM.corrected std.err 
%[1,]         0.0357 0.0282       0.0652   0.240
%[2,]         0.1346 0.1519       0.1819   0.157
%[3,]         0.2221 0.1795       0.1997   0.147
%[4,]         0.3048 0.2952       0.2949   0.153
%[5,]         0.3910 0.3668       0.3579   0.157
%[6,]         0.4860 0.5143       0.4819   0.178
%[7,]         0.6277 0.6556       0.6077   0.208
%
%    plot(cal)
%

The estimated calibration curves are shown in Figure \ref{cal.cox}.
Bias--corrected calibration is very good except for the two groups with
extremely bad prognosis --- their survival is slightly better than
predicted, consistent with regression to the mean.  Even there, the
absolute error is low despite a large relative error.  Hence for this
example it may not be worthwhile to develop a model using shrinkage.
\fig{cal.cox}{Bootstrap estimate of calibration accuracy for 5--year
estimates from the final Cox model.  Dots correspond to apparent
predictive accuracy.  {\bf\tt X} marks the bootstrap--corrected
estimates.}

Now compare this analysis with three previous analyses of this
dataset.  In all three analyses, all continuous covariables were
arbitrarily categorized into intervals and scored with somewhat
arbitrary category codes.  In none of the three were \Co{sbp, dbp,
ekg, ap, bm} \ considered.  Patients having missing values on any of the
candidate predictors were excluded from consideration.

Turn first to Byar and Green \cite{bya80cho}, who used an exponential
survival model and dichotomized treatment by combining placebo and low
dose and combining the two highest doses.  The important predictors
were found to be \Co{hx, sg, sz, hg}, and the following interactions
were detected in an exploratory analysis which did not control for
multiple comparisons: \Co{rx} $\times$ \Co{sg} and \co{rx} $\times\ $
\Co{age}.  These interactions were not significant in the present
model (even if dose were re--coded as in Byar and Green).

Kay \cite{kay86tre} considered Cox models for various causes of death.
For time until all--cause mortality, Kay found that the most important
predictors were \Co{sz, hx, sg, age}. The treatment along with
\Co{age, hx} were significant predictors of cardiovascular death.  The
treatment (in the opposite direction), and \Co{hg,sz,sg} predicted
cancer death.  Treatment and \Co{age, wt} predicted time until death
from other causes.

Sauerbrei and Schumacher \cite{sau92boo} also used a Cox model and
an approach in which a backward elimination procedure was done for
each of 100 bootstrap samples.  The relative frequency of selection of
variables as ``important'' was used as the criterion for inclusion of
variables in the final model.  Variables were retained if they were
selected $\geq 70$ times.  All candidate predictors met this
criterion.  Treatment interactions involving \Co{age} and \Co{sg} were
the most common interactions (56 and 48 bootstrap repetitions,
respectively), but they did not meet the criterion for selection.  The
authors noted that these interactions were misleadingly more
significant in a model which only adjusted for ``significant''
predictors instead of all candidate predictors.

None of the three references just cited provided a model validation or
quantified the predictive discrimination of the final model.

\section{Summary}
Methods were described for developing clinical multivariable
prognostic models and for assessing their calibration and
discrimination. A detailed examination of model assumptions and an
unbiased assessment of predictive accuracy will uncover problems that
may make clinical prediction models misleading or invalid. The
modeling strategy presented in Section \ref{strategy} provides one
sequence of steps for avoiding the pitfalls of multivariable modeling
so that its many advantages can be realized.

\bibliographystyle{ieeetr}
\bibliography{/home/harrelfe/bib/feh}
\end{document}











