Generic function for extracting the predictive partial log-likelihood

Description

Generic function for extracting th predictive partial log-likelihood from a fitted survival model.

Usage

PLL(object, newdata, newtime, newstatus, ...)

Arguments

Details

The predictive partial log-likelihood measures the prediction performance of each model fitted in a boostrap sample, using the data not in this sample. Multiplying by (-2) leads to a deviance-like measure, which means that small values indicate good prediction performance.

peperr requires function PLL.class in case of survival response, for each model fit of class class. Available methods include PLL.CoxBoost, PLL.coxph, PLL.penfit, PLL.coxnet, PLL.grpsurv, PLL.ncvsurv, PLL.mboost_coxph, and PLL.SGL_cox.

Value

Numeric predictive partial log-likelihood returned by the method implementation.

Predictive partial log-likelihood for a CoxBoost model fit

Description

Extracts the predictive partial log-likelihood from a CoxBoost model fit.

Usage

## S3 method for class 'CoxBoost'
PLL(object, newdata, newtime, newstatus, complexity, ...)

Arguments

Details

Used by function peperr, if function fit.CoxBoost is used for model fit.

Since CoxBoost is only suggested by peperr, install it before calling this function.

Value

Numeric predictive partial log-likelihood returned by predict.CoxBoost.

Predictive partial log-likelihood for Cox poportional hazards model

Description

Extracts the predictive partial log-likelihood from a coxph model fit.

Usage

## S3 method for class 'coxph'
PLL(object, newdata, newtime, newstatus, complexity, ...)

Arguments

Details

Used by function peperr, if function fit.coxph is used for model fit.

Value

Vector of length n_new

Predictive partial log-likelihood for penalized Cox models

Description

Computes the predictive partial log-likelihood for penfit survival models.

Usage

## S3 method for class 'penfit'
PLL(object, newdata, newtime, newstatus, complexity = NULL, ...)

Arguments

Value

Numeric predictive partial log-likelihood.

See Also

PLL, fit.penalized, predictProb.penfit

Determine the Brier score for a fitted model

Description

Evaluate the Brier score, i.e. prediction error, for a fitted model on new data. To be used as argument aggregation.fun in peperr call.

Usage

aggregation.brier(full.data=NULL, response, x, model, cplx=NULL,  
type=c("apparent", "noinf"), fullsample.attr = NULL, ...)

Arguments

Details

The empirical Brier score is the mean of the squared difference of the risk prediction and the true value of all observations and takes values between 0 and 1, where small values indicate good prediction performance of the risk prediction model.

Value

Scalar, indicating the empirical Brier score.

Determine the missclassification rate for a fitted model

Description

Evaluate the misclassification rate, i.e. prediction error, for a fitted model on new data. To use as argument aggregation.fun in peperr call.

Usage

aggregation.misclass(full.data=NULL, response, x, model, cplx=NULL,  
type=c("apparent", "noinf"), fullsample.attr = NULL, ...)

Arguments

Details

Misclassification rate is the ratio of observations for which prediction of response is wrong.

Value

Scalar, indicating the misclassification rate.

Determine the prediction error curve for a fitted model

Description

Interface to pmpec, for conforming to the structure required by the argument aggregation.fun in peperr call. Evaluates the prediction error curve, i.e. the Brier score tracked over time, for a fitted survival model.

Usage

aggregation.pmpec(full.data, response, x, model, cplx=NULL, times = NULL, 
   type=c("apparent", "noinf"), fullsample.attr = NULL, ...)

Arguments

Details

If no evaluation time points are passed, they are generated using all uncensored time points if their number is smaller than 100, or 100 time points up to the 95% quantile of the uncensored time points are taken.

pmpec requires a predictProb method for the class of the fitted model, i.e. for a model of class class predictProb.class.

Value

A matrix with one row. Each column represents the estimated prediction error of the fit at the time points.

See Also

peperr, predictProb, pmpec

Complexity selection helpers for penalized backends

Description

Wrapper functions for choosing the penalized tuning parameter used by fit.penalized, fit.LASSO, and fit.fusedLASSO.

Usage

complexity.cv.penalized(response, x, full.data, ...)
complexity.LASSO(response, x, full.data, ...)
complexity.fusedLASSO(response, x, full.data, ...)

Arguments

Details

complexity.cv.penalized profiles the cross-validated lambda1 values returned by profL1 and keeps the maximizer. complexity.LASSO is a compatibility alias. complexity.fusedLASSO performs the same search while forwarding fusedl = TRUE.

Value

A scalar tuning value, typically the selected lambda1.

See Also

fit.penalized, fit.LASSO, fit.fusedLASSO, profL1

Interface function for CoxBoost complexity selection via integrated prediction error curves

Description

Determines the number of boosting steps for a survival model fitted by CoxBoost via integrated prediction error curve estimates, conforming to the calling convention required by argument complexity in peperr call.

Usage

complexity.ipec.CoxBoost(response, x, boot.n.c = 10, boost.steps = 100,
   eval.times = NULL, smooth = FALSE, full.data, ...)

complexity.ripec.CoxBoost(response, x, boot.n.c = 10, boost.steps = 100,
   eval.times = NULL, smooth = FALSE, full.data, ...)

Arguments

Details

Plotting the .632+ estimator for each time point given in eval.times results in a prediction error curve. A summary measure can be obtained by integrating over time. complexity.ripec.CoxBoost computes a Riemann integral, while complexity.ipec.CoxBoost uses a Lebesgue-like integral taking Kaplan-Meier estimates as weights. The number of boosting steps for which the minimal integrated error is obtained is returned.

Since CoxBoost is only suggested by peperr, install it before calling these functions. If smooth = TRUE, locfit is also required.

Value

Scalar value giving the number of boosting steps.

See Also

peperr, CoxBoost

Interface for CoxBoost selection of the optimal number of boosting steps via cross-validation

Description

Determines the number of boosting steps for a survival model fitted by CoxBoost via cross-validation, conforming to the calling convention required by argument complexity in peperr call.

Usage

complexity.mincv.CoxBoost(response, x, full.data, ...)

Arguments

Details

Function is basically a wrapper around cv.CoxBoost of package CoxBoost. A K-fold cross-validation (default K = 10) is performed to search the optimal number of boosting steps, per default in the interval [0, maxstepno]. Calling peperr, the default arguments of cv.CoxBoost can be changed by passing a named list containing these as argument args.complexity.

Since CoxBoost is only suggested by peperr, install it before calling this function.

Value

Scalar value giving the optimal number of boosting steps.

See Also

peperr, cv.CoxBoost

Extract functions, libraries and global variables to be loaded onto a compute cluster

Description

Automatic extraction of functions, libraries and global variables employed passed functions. Designed for peperr call, see Details section there.

Usage

extract.fun(funs = NULL)

Arguments

Details

This function is necessary for compute cluster situations where for computation on nodes required functions, libraries and variables have to be loaded explicitly on each node. Avoids loading of whole global environment which might include the unnecessary loading of huge data sets.

It might have problems in some cases, especially it is not able to extract the library of a function that has no namespace. Similarly, it can only extract a required library if it is loaded, or if the function contains a require or library call.

Value

list containing

See Also

peperr

Examples

# 1. Simplified example for illustration
## Not run: 
library(CoxBoost)
a <- function(){
# some calculation
}

b<-function(){
# some other calculation
x <- cv.CoxBoost()
# z is global variable
y <- a(z)
}

# list with packages, functions and variables required for b:
extract.fun(list(b))

# 2. As called by default in peperr example
extract.fun(list(fit.CoxBoost, aggregation.pmpec))

## End(Not run)

Interface function for fitting a CoxBoost model

Description

Interface for fitting survival models by CoxBoost, conforming to the requirements for argument fit.fun in peperr call.

Usage

fit.CoxBoost(response, x, cplx, ...)

Arguments

Details

Function is basically a wrapper around CoxBoost of package CoxBoost. A Cox proportional hazards model is fitted by componentwise likelihood based boosting, especially suited for models with many covariates and few observations.

Since CoxBoost is only suggested by peperr, install it before calling this function.

Value

CoxBoost object

See Also

peperr, CoxBoost

Sparse-group-lasso Cox backend

Description

Wrappers for Cox models fitted with SGL.

Usage

fit.SGL.cox(response, x, cplx, index, ...)
complexity.cvSGL.cox(response, x, full.data, index, ...)
## S3 method for class 'SGL_cox'
predictProb(object, response, x, times, complexity = NULL, ...)
## S3 method for class 'SGL_cox'
PLL(object, newdata, newtime, newstatus, complexity = NULL, ...)

Arguments

Details

The SGL predictor returns relative risks. The peperr wrapper reconstructs survival probabilities from those risks with a Breslow baseline estimated on the stored training data.

Value

Fitted SGL_cox objects, selected lambda values, survival-probability matrices, and numeric predictive partial log-likelihood values, respectively.

See Also

peperr, SGL, cvSGL

Interface function for fitting a Cox proportional hazards model

Description

Interface for fitting survival models by Cox proporional hazards model, conforming to the requirements for argument fit.fun in peperr call.

Usage

fit.coxph(response, x, cplx, ...)

Arguments

Details

Function is basically a wrapper around coxph of package survival.

Value

coxph object

See Also

peperr, coxph

glmnet Cox backend

Description

Wrappers that integrate glmnet Cox models into the peperr fit/complexity/prediction interface.

Usage

fit.glmnet(response, x, cplx, ...)
complexity.cv.glmnet(response, x, full.data, ...)
## S3 method for class 'coxnet'
predictProb(object, response, x, times, complexity = NULL, ...)
## S3 method for class 'coxnet'
PLL(object, newdata, newtime, newstatus, complexity = NULL, ...)

Arguments

Details

The backend stores the training design matrix and survival outcome on the fitted object so that the survfit.coxnet method can be reused later when pmpec calls predictProb.

Value

fit.glmnet returns a fitted coxnet object. complexity.cv.glmnet returns a scalar lambda. predictProb.coxnet returns an n * length(times) survival-probability matrix. PLL.coxnet returns a numeric predictive partial log-likelihood.

See Also

peperr, predictProb, PLL, glmnet, cv.glmnet

grpreg grouped Cox backend

Description

Wrappers for grouped survival models fitted with grpreg.

Usage

fit.grpsurv(response, x, cplx, group, ...)
complexity.cv.grpsurv(response, x, full.data, group, ...)
## S3 method for class 'grpsurv'
predictProb(object, response, x, times, complexity = NULL, ...)
## S3 method for class 'grpsurv'
PLL(object, newdata, newtime, newstatus, complexity = NULL, ...)

Arguments

Details

The grouped backend exposes grouped lasso, grouped MCP, grouped SCAD, and related penalties through the native grpreg arguments passed in ....

Value

Fitted grpsurv objects, scalar lambda values, survival-probability matrices, and numeric predictive partial log-likelihood values, respectively.

See Also

peperr, grpsurv, cv.grpreg

mboost Cox backend

Description

Wrappers for Cox boosting models fitted with mboost.

Usage

fit.mboost.coxph(response, x, cplx, ...)
complexity.cvrisk.mboost(response, x, full.data, max.mstop = 100L, folds = NULL, ...)
## S3 method for class 'mboost_coxph'
predictProb(object, response, x, times, complexity = NULL, ...)
## S3 method for class 'mboost_coxph'
PLL(object, newdata, newtime, newstatus, complexity = NULL, ...)

Arguments

Details

If mboost provides native survival-curve predictions for the fitted model, predictProb.mboost_coxph uses them. Otherwise it falls back to a Breslow baseline estimated from the stored training data and linear predictors.

Value

Fitted mboost_coxph objects, selected stopping iterations, survival-probability matrices, and numeric predictive partial log-likelihood values, respectively.

See Also

peperr, glmboost, cvrisk

ncvreg survival backend

Description

Wrappers for Cox models fitted with ncvreg.

Usage

fit.ncvsurv(response, x, cplx, ...)
complexity.cv.ncvsurv(response, x, full.data, ...)
## S3 method for class 'ncvsurv'
predictProb(object, response, x, times, complexity = NULL, ...)
## S3 method for class 'ncvsurv'
PLL(object, newdata, newtime, newstatus, complexity = NULL, ...)

Arguments

Details

The implementation first tries the native ncvreg survival prediction method. If that is unavailable for the installed version, it falls back to a Breslow baseline built from stored training data and linear predictors.

Value

Fitted ncvsurv objects, selected lambda values, survival-probability matrices, and numeric predictive partial log-likelihood values, respectively.

See Also

peperr, predictProb, PLL

Penalized backend fitting wrappers

Description

Interface functions for fitting penalized models in peperr. The generic wrapper is fit.penalized; fit.LASSO and fit.fusedLASSO are convenience aliases for plain and fused lasso fits.

Usage

fit.penalized(response, x, cplx, ...)
fit.LASSO(response, x, cplx, ...)
fit.fusedLASSO(response, x, cplx, ...)

Arguments

Details

fit.penalized is the general wrapper around penalized. For survival models it also stores the training design matrix and response so that predictProb.penfit and PLL.penfit can be used later by pmpec and peperr.

fit.LASSO simply calls fit.penalized. fit.fusedLASSO calls the same backend with fusedl = TRUE.

Value

A fitted penfit object.

See Also

complexity.cv.penalized, predictProb.penfit, PLL.penfit, penalized

randomForestSRC survival backend

Description

Wrappers for random survival forests fitted with randomForestSRC.

Usage

fit.rfsrc(response, x, cplx, ...)
complexity.oob.rfsrc(response, x, full.data,
    mtry = unique(pmax(1L, c(floor(sqrt(ncol(x))), floor(ncol(x)/3), ncol(x)))),
    nodesize = c(5L, 15L), ntree = 200L, ...)
## S3 method for class 'rfsrc'
predictProb(object, response, x, times, complexity = NULL, ...)

Arguments

Details

No PLL.rfsrc method is provided, because random survival forests do not naturally expose a Cox-style partial likelihood. They can still be used with pmpec through predictProb.rfsrc.

Value

Fitted rfsrc objects, tuning lists containing mtry and nodesize, and survival-probability matrices.

See Also

peperr, rfsrc, predict.rfsrc

Integrated prediction error curve

Description

Summary measures of prediction error curves

Usage

ipec(pe, eval.times, type=c("Riemann", "Lebesgue", "relativeLebesgue"), response=NULL)

Arguments

Details

For survival data, prediction error at each evaluation time point can be extracted of a peperr object by function perr. A summary measure can then be obtained via intgrating over time. Note that the time points used for evaluation are stored in list element attribute of the peperr object.

Value

See Also

perr

Examples

## Not run: 
n <- 200
p <- 100
beta <- c(rep(1,10),rep(0,p-10))
x <- matrix(rnorm(n*p),n,p)
real.time <- -(log(runif(n)))/(10*exp(drop(x %*% beta)))
cens.time <- rexp(n,rate=1/10)
status <- ifelse(real.time <= cens.time,1,0)
time <- ifelse(real.time <= cens.time,real.time,cens.time)

# Example:
# Obtain prediction error estimate fitting a Cox proportional hazards model
# using CoxBoost 
# through 10 bootstrap samples 
# with fixed complexity 50 and 75
# and aggregate using prediction error curves
peperr.object <- peperr(response=Surv(time, status), x=x, 
   fit.fun=fit.CoxBoost, complexity=c(50, 75), 
   indices=resample.indices(n=length(time), method="sub632", sample.n=10))
# 632+ estimate for both complexity values at each time point
prederr <- perr(peperr.object)
# Integrated prediction error curve for both complexity values
ipec(prederr, eval.times=peperr.object$attribute, response=Surv(time, status))

## End(Not run)

Parallelised Estimation of Prediction Error

Description

Prediction error estimation for regression models via resampling techniques. Potentially parallelised, if compute cluster is available.

Usage

peperr(response, x, 
	indices = NULL, 
	fit.fun, complexity = NULL, args.fit = NULL, args.complexity = NULL,
	parallel = NULL, cpus = 2, clustertype=NULL, clusterhosts=NULL,
	noclusterstart = FALSE, noclusterstop=FALSE,
	aggregation.fun=NULL, args.aggregation = NULL, 
	load.list = extract.fun(list(fit.fun, complexity, aggregation.fun)),
	load.vars = NULL, load.all = FALSE, 
	trace = FALSE, debug = FALSE,
	peperr.lib.loc=NULL, 
        RNG=c("RNGstream", "SPRNG", "fixed", "none"), seed=NULL, 
        lb=FALSE, sr=FALSE, sr.name="default", sr.restore=FALSE)

Arguments

Details

Validation of new model fitting approaches requires the proper use of resampling techniques for prediction error estimation. Especially in high-dimensional data situations the computational demand might be huge. peperr accelerates computation through automatically parallelisation of the resampling procedure, if a compute cluster is available. A noticeable speed-up is reached even when using a dual-core processor.

Resampling based prediction error estimation requires for each split in training and test data the following steps: a) selection of model complexity (if desired), using the training data set, b) fitting the model with the selected (or a given) complexity on the training set and c) measurement of prediction error on the corresponding test set.

Functions for fitting the model, determination of model complexity, if required by the fitting procedure, and aggregating the prediction error are passed as arguments fit.fun, complexity and aggregation.fun. Already available functions are

for model fit: fit.CoxBoost, fit.coxph, fit.penalized, fit.LASSO, fit.fusedLASSO, fit.glmnet, fit.grpsurv, fit.ncvsurv, fit.rfsrc, fit.mboost.coxph, fit.SGL.cox

to determine complexity: complexity.mincv.CoxBoost, complexity.ipec.CoxBoost, complexity.ripec.CoxBoost, complexity.cv.penalized, complexity.LASSO, complexity.fusedLASSO, complexity.cv.glmnet, complexity.cv.grpsurv, complexity.cv.ncvsurv, complexity.oob.rfsrc, complexity.cvrisk.mboost, complexity.cvSGL.cox

to aggregate prediction error: aggregation.pmpec, aggregation.brier, aggregation.misclass

Function peperr is especially designed for evaluation of newly developed model fitting routines. For that, own routines can be passed as arguments to the peperr call. They are incorporated as follows (also compare existing functions, as named above):

1. Model fitting techniques, which require selection of one or more complexity parameters, often provide routines based on cross-validation or similar to determine this parameter. If this routine is already at hand, the complexity function needed for the peperr call is not more than a wrapper around that, which consists of providing the data in the required way, calling the routine and return the selected complexity value(s).

2. For a given model fitting routine the fitting function, which is passed to the peperr call as argument fit.fun, is not more than a wrapper around that. Explicitly, response and matrix of covariates have to be transformed to the required form, if necessary, the routine is called with the passed complexity value, if required, and the fitted prediction model is returned.

3. Prediction error is estimated using a fitted model and a data set, by any kind of comparison of the true and the predicted response values. In case of survival response, apparent error (type apparent), which means that the prediction error is estimated in the same data set as used for model fitting, and no-information error (type noinf), which calculates the prediction error in permuted data, have to be provided. Note that the aggregation function returns the error with an additional attribute called addattr. The evaluation time points have to be stored there to allow later access.

4. In case of survival response, the user may additionally provide a function for partial log likelihood calculation, if he uses an own function for model fit, called PLL.class. If prediction error curves are used for aggregation (aggregation.pmpec), a predictProb method has to be provided, i.e. for each model of class class predictProb.class, see there.

Concerning parallelisation, there are three possibilities to run peperr:

• Start R on commandline with sfCluster and preferred options, for example number of cpus. Leave the three arguments parallel, clustertype and nodes unchanged.

• Use any other cluster solution supported by snowfall, i.e. LAM/MPI, socket, PVM, NWS (set argument clustertype). Argument parallel has to be set to TRUE and number of cpus can be chosen by argument nodes)

• If no cluster is used, R works sequentially. Keep parallel=NULL. No parallelisation takes place and therefore no speed up can be obtained.

In general, if parallel=NULL, all information concerning the cluster set-up is taken from commandline, else, it can be specified using the three arguments parallel, clustertype, nodes, and, if necessary, clusterhosts.

sfCluster is a Unix tool for flexible and comfortable managment of parallel R processes. However, peperr is usable with any other cluster solution supported by snowfall, i.e. sfCluster has not to be installed to use package peperr. Note that this may require cluster handling by the user, e.g. manually shut down with 'lamhalt' on commandline for type="MPI". But, using a socket cluster (argument parallel=TRUE and clustertype="SOCK"), does not require any extra installation.

Note that the run time cannot speed up anymore if the number of nodes is chosen higher than the number of passed training/test samples plus one, as parallelisation takes place in the resampling procedure and one additional run is used for computation on the full sample.

If not running in sequential mode, a specified number of R processes called nodes is spawned for parallel execution of the resampling procedure (see above). This requires to provide all variables, functions and libraries necessary for computation on each of these R processes, so explicitly all variables, functions and libraries required by the, potentially user-defined, functions fit.fun, complexity and aggregation.fun. The simplest possibility is to load the whole content of the global environment on each node and all loaded libraries. This is done by setting argument load.all=TRUE. This is not the default, as a huge amount of data is potentially loaded to each node unnecessarily. Function extract.fun is provided to extract the functions and libraries needed, automatically called at each call of function peperr. Note that all required libraries have to be located in the standard library search path (obtained by .libPaths()). Another alternative is to load required data manually on the slaves, using snowfall functions sfLibrary, sfExport and sfExportAll. Then, argument noclusterstart has to be switched to TRUE. Additionally, argument load.list could be set to NULL, to avoid potentially overwriting of functions and variables loaded to the cluster nodes automatically.

Note that a set.seed call before calling function peperr is not sufficient to allow reproducibility of results when running in parallel mode, as the slave R processes are not affected as they are own R instances. peperr provides two possibilities to make results reproducible:

• Use RNG="RNGstream" for independent parallel random number streams initialized on the cluster nodes via sfClusterSetupRNG from package snowfall. This requires package rlecuyer. A seed can be specified using argument seed, else the default values are taken. A set.seed call on the master is required additionally and argument lb=FALSE, see below.

• RNG="SPRNG" is kept as a recognized legacy value so existing code fails with an informative message, but it is no longer supported because rsprng has been removed from CRAN.

• If RNG="fixed", a seed has to be specified. This can be either an integer or a vector of length number of samples +2. In the second case, the first entry is used for the main R process, the next number of samples ones for each sample run (in parallel execution mode on slave R processes) and the last one for computation on full sample (as well on slave R process in parallel execution mode). Passing integer x is equivalent to passing vector x+(0:(number of samples+1)). This procedure allows reproducibility in any case, i.e. also if the number of parallel processes changes as well as in sequential execution.

Load balancing (argument lb) means, that a slave gets a new job immediately after the previous is finished. This speeds up computation, but may change the order of jobs. Due to that, results are only reproducible, if RNG="fixed" is used.

Value

Object of class peperr

Author(s)

Christine Porzelius cp@fdm.uni-freiburg.de, Harald Binder

References

Binder, H. and Schumacher, M. (2008) Adapting prediction error estimates for biased complexity selection in high-dimensional bootstrap samples. Statistical Applications in Genetics and Molecular Biology, 7:1.

Porzelius, C., Binder, H., Schumacher, M. (2008) Parallelised prediction error estimation for evaluation of high-dimensional models. Manuscript.

See Also

perr, resample.indices, extract.fun

Examples

# Generate survival data with 10 informative covariates
n <- 200
p <- 100
beta <- c(rep(1,10),rep(0,p-10))
x <- matrix(rnorm(n*p),n,p)
real.time <- -(log(runif(n)))/(10*exp(drop(x %*% beta)))
cens.time <- rexp(n,rate=1/10)
status <- ifelse(real.time <= cens.time,1,0)
time <- ifelse(real.time <= cens.time,real.time,cens.time)

## Not run: 
# A: R runs sequential or R is started on commandline with desired options 
# (for example using sfCluster: sfCluster -i --cpus=5)
# Example A1:
# Obtain prediction error estimate fitting a Cox proportional hazards model
# using CoxBoost 
# through 10 bootstrap samples 
# with fixed complexity 50 and 75
# and aggregate using prediction error curves (default setting)

peperr.object1 <- peperr(response=Surv(time, status), x=x, 
   fit.fun=fit.CoxBoost, complexity=c(50, 75), 
   indices=resample.indices(n=length(time), method="sub632", sample.n=10))
peperr.object1

# Diagnostic plots
plot(peperr.object1)

# Extraction of prediction error curves (.632+ prediction error estimate), 
# blue line corresponds to complexity 50, 
# red one to complexity 75
plot(peperr.object1$attribute,
   perr(peperr.object1)[1,], type="l", col="blue",
   xlab="Evaluation time points", ylab="Prediction error")
lines(peperr.object1$attribute, 
   perr(peperr.object1)[2,], col="red")

# Example A2:
# As Example A1, but
# with complexity selected through a cross-validation procedure
# and extra argument 'penalty' passed to fit function and complexity function
peperr.object2 <- peperr(response=Surv(time, status), x=x, 
   fit.fun=fit.CoxBoost, args.fit=list(penalty=100),
   complexity=complexity.mincv.CoxBoost, args.complexity=list(penalty=100),
   indices=resample.indices(n=length(time), method="sub632", sample.n=10),
   trace=TRUE)
peperr.object2

# Diagnostic plots
plot(peperr.object2)

# Example A3:
# As Example A2, but
# with extra argument 'times', specifying the evaluation times passed to aggregation.fun
# and seed, for reproducibility of results
# Note: set.seed() is required additional to argument 'seed', 
# as function 'resample.indices' is used in peperr call.
set.seed(123)
peperr.object3 <- peperr(response=Surv(time, status), x=x, 
   fit.fun=fit.CoxBoost, args.fit=list(penalty=100),
   complexity=complexity.mincv.CoxBoost, args.complexity=list(penalty=100),
   indices=resample.indices(n=length(time), method="sub632", sample.n=10),
   args.aggregation=list(times=seq(0, quantile(time, probs=0.9), length.out=100)),
   trace=TRUE, RNG="fixed", seed=321)
peperr.object3

# Diagnostic plots
plot(peperr.object3)

# B: R is started sequential, desired cluster options are given as arguments
# Example B1:
# As example A1, but using a socket cluster and 3 CPUs
peperr.object4 <- peperr(response=Surv(time, status), x=x, 
   fit.fun=fit.CoxBoost, complexity=c(50, 75), 
   indices=resample.indices(n=length(time), method="sub632", sample.n=10),
   parallel=TRUE, clustertype="SOCK", cpus=3)

## End(Not run)

Prediction error estimates

Description

Extracts prediction error estimates from peperr objects.

Usage

perr(peperrobject, 
    type = c("632p", "632", "apparent", "NoInf", "resample", "nullmodel"))

Arguments

Details

The .632 and the .632+ prediction error estimates are weighted combinations of the apparent error and bootstrap cross-validation error estimate, for survival data at given time points.

Value

If type="632p" or type="632": Prediction error: Matrix, with one row per complexity value.

If type="apparent": Apparent error of the full data set. Matrix: One row per complexity value. In case of survival response, columns correspond to evaluation timepoints, which are given in attribute addattr.

If type="NoInf": No-information error of the full data set, i. e. evaluation in permuted data. Matrix: One row per complexity value. Columns correspond to evaluation timepoints, which are given in attribute addattr.

If type="resample": Matrix. Mean prediction error of resampling test samples, one row per complexity value.

If type="nullmodel": Vector or scalar: Null model prediction error, i.e. of fit without information of covariates. In case of survival response Kaplan-Meier estimate at each time point, if response is binary logistic regression model, else not available.

References

Binder, H. and Schumacher, M. (2008) Adapting prediction error estimates for biased complexity selection in high-dimensional bootstrap samples. Statistical Applications in Genetics and Molecular Biology, 7:1.

Gerds, T. and Schumacher, M. (2007) Efron-type measures of prediction error for survival analysis. Biometrics, 63, 1283–1287.

Schumacher, M. and Binder, H., and Gerds, T. (2007) Assessment of Survival Prediction Models in High-Dimensional Settings. Bioinformatics, 23, 1768-1774.

See Also

peperr, ipec

Examples

## Not run: 
n <- 200
p <- 100
beta <- c(rep(1,10),rep(0,p-10))
x <- matrix(rnorm(n*p),n,p)
real.time <- -(log(runif(n)))/(10*exp(drop(x %*% beta)))
cens.time <- rexp(n,rate=1/10)
status <- ifelse(real.time <= cens.time,1,0)
time <- ifelse(real.time <= cens.time,real.time,cens.time)

# Example:
# Obtain prediction error estimate fitting a Cox proportional hazards model
# using CoxBoost 
# through 10 bootstrap samples 
# with fixed complexity 50 and 75
# and aggregate using prediction error curves
peperr.object <- peperr(response=Surv(time, status), x=x, 
   fit.fun=fit.CoxBoost, complexity=c(50, 75), 
   indices=resample.indices(n=length(time), method="sub632", sample.n=10))
# 632+ estimate for both complexity values at each time point
perr(peperr.object)


## End(Not run)

Plot method for peperr object

Description

Plots, allowing to get a first impression of the prediction error estimates and to check complexity selection in bootstrap samples.

Usage

## S3 method for class 'peperr'
plot(x, y, ...)

Arguments

Details

The plots provide a simple and fast overview of the results of the estimation of the prediction error through resampling. Which plots are shown depends on if complexity was selected, i.e., a function was passed in the peperr call for complexity, or explicitly passed. In case of survival response, prediction error curves are shown. In case of binary response, where one complexity value is passed explicitly, no plot is available. Especially in the case that complexity is selected in each bootstrap sample, these diagnostic plots help to check whether the resampling procedure works adequately and to detect specific problems due to high-dimensional data structures.

Examples

## Not run: 
n <- 200
p <- 100
beta <- c(rep(1,10),rep(0,p-10))
x <- matrix(rnorm(n*p),n,p)
real.time <- -(log(runif(n)))/(10*exp(drop(x %*% beta)))
cens.time <- rexp(n,rate=1/10)
status <- ifelse(real.time <= cens.time,1,0)
time <- ifelse(real.time <= cens.time,real.time,cens.time)

peperr.object1 <- peperr(response=Surv(time, status), x=x, 
   fit.fun=fit.CoxBoost, complexity=c(50, 75), 
   indices=resample.indices(n=length(time), method="sub632", sample.n=10))
plot(peperr.object1)

peperr.object2 <- peperr(response=Surv(time, status), x=x, 
   fit.fun=fit.CoxBoost, args.fit=list(penalty=100),
   complexity=complexity.mincv.CoxBoost, args.complexity=list(penalty=100),
   indices=resample.indices(n=length(time), method="sub632", sample.n=10),
   trace=TRUE)
plot(peperr.object2)

peperr.object3 <- peperr(response=Surv(time, status), x=x, 
   fit.fun=fit.CoxBoost, args.fit=list(penalty=100),
   complexity=complexity.mincv.CoxBoost, args.complexity=list(penalty=100),
   indices=resample.indices(n=length(time), method="sub632", sample.n=10),
   args.aggregation=list(times=seq(0, quantile(time, probs=0.9), length.out=100)),
   trace=TRUE)
plot(peperr.object3)

## End(Not run)

Calculate prediction error curves

Description

Calculation of prediction error curve from a survival response and predicted probabilities of survival.

Usage

pmpec(object, response=NULL, x=NULL, times, model.args=NULL, 
    type=c("PErr","NoInf"), external.time=NULL, external.status=NULL, 
    data=NULL)

Arguments

Details

Prediction error of survival data is measured by the Brier score, which considers the squared difference of the true event status at a given time point and the predicted event status by a risk prediction model at that time. A prediction error curve is the weighted mean Brier score as a function of time at time points in times (see References).

pmpec requires a predictProb method for the class of the fitted model, i.e. for a model of class class predictProb.class.

pmpec is implemented to behave similar to function pec of package pec, which provides several predictProb methods.

In bootstrap framework, data contains only a part of the full data set. For censoring distribution, the full data should be used to avoid extreme variance in case of small data sets. For that, the observed times and status values can be passed as argument external.time and external.status.

Value

Vector of prediction error estimates at each time point given in time.

Author(s)

Harald Binder

References

Gerds, A. and Schumacher, M. (2006) Consistent estimation of the expected Brier score in general survival models with right-censored event times. Biometrical Journal, 48, 1029–1040.

Schoop, R. (2008) Predictive accuracy of failure time models with longitudinal covariates. PhD thesis, University of Freiburg. http://www.freidok.uni-freiburg.de/volltexte/4995/.

See Also

predictProb, pec

Generic function for extracting predicted survival probabilities

Description

Generic function for extraction of predicted survival probabilities from a fitted survival model conforming to the interface required by pmpec.

Usage

predictProb(object, response, x, ...)

Arguments

Details

pmpec requires a predictProb.class function for each model fit of class class. It extracts the predicted probability of survival from this model.

See existing predictProb functions, at the time predictProb.CoxBoost, predictProb.coxph, predictProb.survfit, predictProb.penfit, predictProb.coxnet, predictProb.grpsurv, predictProb.ncvsurv, predictProb.rfsrc, predictProb.mboost_coxph, and predictProb.SGL_cox.

If desired predictProb function for class class is not available in peperr, but implemented in package pec as predictSurvProb.class, it can easily be transformed as predictProb method.

Value

Matrix with predicted probabilities for each evaluation time point in times (columns) and each new observation (rows).

Extract predicted survival probabilities from a CoxBoost fit

Description

Extracts predicted survival probabilities from a survival model fitted by CoxBoost, providing an interface as required by pmpec.

Usage

## S3 method for class 'CoxBoost'
predictProb(object, response, x, times, complexity, ...)

Arguments

Details

Since CoxBoost is only suggested by peperr, install it before calling this function.

Value

Matrix with probabilities for each evaluation time point in times (columns) and each new observation (rows).

Extract predicted survival probabilities from a coxph object

Description

Extracts predicted survival probabilities for survival models fitted by Cox proportional hazards model, providing an interface as required by pmpec.

Usage

## S3 method for class 'coxph'
predictProb(object, response, x, times, ...)

Arguments

Value

Matrix with probabilities for each evaluation time point in times(columns) and each new observation (rows).

Predicted survival probabilities for penalized Cox models

Description

Extracts survival probabilities from a penfit survival model for use in pmpec.

Usage

## S3 method for class 'penfit'
predictProb(object, response, x, times, complexity = NULL, ...)

Arguments

Value

Matrix of survival probabilities with one row per observation and one column per evaluation time.

See Also

predictProb, fit.penalized, PLL.penfit

Extract predicted survival probabilities from a survfit object

Description

Extracts predicted survival probabilities for survival models fitted by survfit, providing an interface as required by pmpec.

Usage

## S3 method for class 'survfit'
predictProb(object, response, x, times, train.data, ...)

Arguments

Value

Matrix with probabilities for each evaluation time point in times(columns) and each new observation (rows).

Generation of indices for resampling Procedure

Description

Generates training and test set indices for use in resampling estimation of prediction error, e.g. cross-validation or bootstrap (with and without replacement).

Usage

resample.indices(n, sample.n = 100, method = c("no", "cv" ,"boot", "sub632"))

Arguments

Details

As each bootstrap sample should be taken as if new data, complexity selection should be carried out in each bootstrap sample. Binder and Schumacher show that when bootstrap samples are drawn with replacement, often too complex models are obtained in high-dimensional data settings. They recommend to draw bootstrap samples without replacement, each of size round(0.632 * n), which equals the expected number of unique observations in one bootstrap sample drawn with replacement, to avoid biased complexity selection and improve predictive power of the resulting models.

Value

A list containing two lists of length sample.n:

References

Binder, H. and Schumacher, M. (2008) Adapting prediction error estimates for biased complexity selection in high-dimensional bootstrap samples. Statistical Applications in Genetics and Molecular Biology, 7:1.

See Also

peperr

Examples

# generate dataset: 100 patients, 20 covariates
data <- matrix(rnorm(2000), nrow=100)

# generate indices for training and test data for 10-fold cross-validation
indices <- resample.indices(n=100, sample.n = 10, method = "cv")

# create training and test data via indices
trainingsample1 <- data[indices$sample.index[[1]],]
testsample1 <- data[indices$not.in.sample[[1]],]