pec (version 1.0.7)

pec: Prediction error curves

Description

Evaluating the performance of risk prediction models in survival analysis. The Brier score is a weighted average of the squared distances between the observed survival status and the predicted survival probability of a model. Roughly the weights correspond to the probabilities of not being censored and they might depend on covariates. Prediction error curves are obtained when the Brier score is followed over time. Bootstrap-crossvalidation techniques are implemented to estimate the generalization error.

Usage

pec(object,...)
## S3 method for class 'list':
pec(object,
    formula,
    data,
    times,
    start,
    maxtime,
    exact=TRUE,
    exactness=100,
    fillChar=NA,
    cens.model="cox",
    replan="none",
    B,
    M,
    model.args=NULL,
    model.parms=NULL,
    keep.matrix=FALSE,
    import=NULL,
    export=NULL,
    na.accept=0,
    verbose=TRUE,
    ...)
## S3 method for class 'coxph':
pec(object,...)
## S3 method for class 'glm':
pec(object,...)
## S3 method for class 'cph':
pec(object,...)
## S3 method for class 'prodlim':
pec(object,...)
## S3 method for class 'survfit':
pec(object,...)
## S3 method for class 'aalen':
pec(object,...)

Arguments

object

A named list of prediction models. Each entry is either an R-object for which a predictSurvProb method exists (see details) or a call that can be evaluated to such an R-object. For resam

formula

A formula. If missing the formula of the first model in list is used. The left hand side is used to finde the status response variable in data. For right censored data, the right hand side of the formula is used to specif

data

A data frame in which to validate the prediction models and to fit the censoring model. If missing the data of the first model in list is used.

times

A vector of timepoints. At each timepoint the prediction error curves are estimated. If missing and exact==TRUE use all the unique values of the response variable. If missing and exact==FALSE use a equidistan

start

Minimal time for estimating the prediction error curves. If missing and formula defines a Surv or Hist object then start defaults to 0, otherwise to the smallest observed value

maxtime

Maximal time for estimating the prediction error curves. If missing the largest value of the response variable is used.

exact

Logical. If TRUE estimate the prediction error curves at all the unique values of the response variable. If times are given and exact=TRUE then the times are merged with the unique values of

exactness

An integer that determines how many equidistant gridpoints are used between start and maxtime. The default is 100.

fillChar

Symbol used to fill-in places where the values of the prediction error curves are not available. The default is NA.

cens.model

Method for estimating inverse probability of censoring weigths: cox: A semi-parametric Cox proportional hazard model is fitted to the censoring times marginal: The Kaplan-Meier estimator for the censorin

replan

Method for estimating the prediction error curves. none:{Assess the models in the same data where they are fitted. }

boot: { Bootstrap the prediction error curves B times WITHOUT cross-

B

Number of bootstrap samples. The default depends on argument replan. When replan in c("OutOfBag","Boot632","Boot632plus" the default is 100. For replan="cvK" B is the number of cross-validation c

M

The size of the bootstrap samples for resampling without replacement.

model.args

Experimental. List of extra arguments that can be passed to the predictSurvProb methods. The list must have an entry for each entry in object.

model.parms

Experimental. List of with exactly one entry for each entry in object. Each entry names parts of the value of the fitted models that should be extracted and added to the value.

keep.matrix

Logical. If TRUE add all B prediction error curves to the output.

import

Experimental. Import the bootstrap matrix from an external file. import is a list of arguments that passed to internal function pecImportIndex

export

Experimental. Export the bootstrap matrix to an external file. export is a list of arguments that passed to internal function pecExportIndex together with the matrix.

na.accept

Only for "OutOfBag" error. The maximal number of failures during training the models to the bootstrap samples. Usually a small number relative to B.

verbose

if TRUE the procedure is reporting details of the progress, e.g. it prints the current step in resampling procedures.

...

Difficult to explain

Value

• A pec object for which print, summary, and plot methods defined. The object includes the following components:
• PredErr
  The estimated prediction error according to the replan. A matrix where each column represents the estimated prediction error of a fit at the time points in time.
• AppErr
  The training error or apparent error obtained when the model(s) are evaluated in the same data where they were trained. Only if replan is one of "NoInf", "cvK", "OutOfBag", "Boot632" or "Boot632plus".
• OutOfBagErr
  The prediction error when the model(s) are trained in the bootstrap sample and evaluated in the data that are not in the bootstrap sample. Only if replan is one of "Boot632" or "Boot632plus". When replan="OutOfBag" then the OutOfBagErr is stored in the component PredErr.
• NoInfErr
  The prediction error when the model(s) are evaluated in the permuted data. Only if replan is one of "OutOfBag", "Boot632", or "Boot632plus". For replan="NoInf" the NoInfErr is stored in the component PredErr.
• weight
  The weight used to linear combine the AppErr and the OutOfBagErr Only if replan is one of "Boot632", or "Boot632plus".
• overfit
  Estimated overfit of the model(s). See Efron & Tibshirani (1997, Journal of the American Statistical Association) and Gerds & Schumacher (2007, Biometrics). Only if replan is one of "Boot632", or "Boot632plus".
• call
  The call that produced the object
• time
  The time points at which the prediction error curves change.
• ipcw.fit
  The fitted censoring model that was used for re-weigthing the Brier score residuals. See Gerds & Schumacher (2006, Biometrical Journal)
• n.risk
  The number of subjects at risk for all time points.
• models
  The prediction models fitted in their own data.
• cens.model
  The censoring models.
• maxtime
  The latest timepoint where the prediction error curves are estimated.
• start
  The earliest timepoint where the prediction error curves are estimated.
• exact
  TRUE if the prediction error curves are estimated at all unique values of the response in the full data.
• method
  The method used for estimation of the prediction error.

Details

Missing data in the response or in the input matrix cause a failure. The status of the continuous response variable at cutpoints (times), ie status=1 if the response value exceeds the cutpoint and status=0 otherwise, is compared to predicted event status probabilities which are provided by the prediction models on the basis of covariates. The comparison is done with the Brier score: the quadratic difference between 0-1 response status and predicted probability. There are two different sources for bias when estimating prediction error in right censored survival problems: censoring and high flexibility of the prediction model. The first is controlled by inverse probability of censoring weighting, the second can be controlled by special Monte Carlo simulation. In each step, the resampling procedures reevaluate the prediction model. Technically this is done by replacing the argument object$call$data by the current subset or bootstrap sample of the full data.

For each prediction model there must be a predictSurvProb method: for example, to assess a prediction model which evaluates to a myclass object one defines a function called predictSurvProb.myclass with arguments object,newdata,cutpoints,train.data,...

Such a function takes the object which was fitted with train.data and derives a matrix with predicted event status probabilities for each subject in newdata (rows) and each cutpoint (column) of the response variable that defines an event status.

Currently, predictSurvProb methods are available for the following R-objects: [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

References

E. Graf et al. (1999), Assessment and comparison of prognostic classification schemes for survival data. Statistics in Medicine, vol 18, pp= 2529--2545.

Efron, Tibshirani (1997) Journal of the American Statistical Association 92, 548--560 Improvement On Cross-Validation: The .632+ Bootstrap Method. Gerds, Schumacher (2006), Consistent estimation of the expected Brier score in general survival models with right-censored event times. Biometrical Journal, vol 48, 1029--1040. Thomas A. Gerds, Martin Schumacher (2007) Efron-Type Measures of Prediction Error for Survival Analysis Biometrics (OnlineEarly Articles). doi:10.1111/j.1541-0420.2007.00832.x

Martin Schumacher, Harald Binder, and Thomas Gerds. Assessment of survival prediction models based on microarray data. Bioinformatics, 23(14):1768-74, 2007.

See Also

plot.pec, summary.pec, R2, crps

Examples

# simulate an artificial data frame
# with survival response and two predictors

set.seed(130971)
N <- 300
X1 <- rnorm(N,mean=10,sd=5)
X2 <- rbinom(N,1,.4)
linPred <- .00001+0.2*X1+2.3*X2
T <- sapply(linPred,function(lp){rexp(n=1,exp(-lp))})
C <- rexp(n=300,exp(-mean(linPred)))
dat <- data.frame(time=pmin(T,C),status=as.numeric(T<=C),X1=X1,X2=X2)

# fit some candidate Cox models and compute the Kaplan-Meier estimate 

Models <- list("Kaplan.Meier"=survfit(Surv(time,status)~1,data=dat),
               "Cox.X1"=coxph(Surv(time,status)~X1,data=dat),
               "Cox.X2"=coxph(Surv(time,status)~X2,data=dat),
               "Cox.X1.X2"=coxph(Surv(time,status)~X1+X2,data=dat))

# compute the apparent prediction error 
PredError <- pec(object=Models,
                  formula=Surv(time,status)~X1+X2,
                  data=dat,
                  exact=TRUE,
                  cens.model="marginal",
                  replan="none",
                  B=0,
                  verbose=TRUE)

print(PredError,times=seq(5,30,5))
summary(PredError)
plot(PredError,xlim=c(0,30))

# compute the .632+ estimate of the generalization error 
set.seed(17100)
PredError.632plus <- pec(object=Models,
                  formula=Surv(time,status)~X1+X2,
                  data=dat,
                  exact=TRUE,
                  cens.model="marginal",
                  replan="boot632plus",
                  B=100,
                  verbose=TRUE)

print(PredError.632plus,times=seq(5,30,5))
summary(PredError.632plus)
plot(PredError.632plus,xlim=c(0,30))