| Title: | Performance Assessment for Continuous and Survival Outcomes |
|---|---|
| Description: | Evaluates prediction accuracy for continuous outcomes and the predictive performance of survival models in a unified framework. Implements the R-squared, L-squared, pseudo R-squared and related explained variation measures of Li and Wang (2019) <doi:10.1080/01621459.2018.1515079>, the concordance indices of Harrell et al. (1982) <doi:10.1001/jama.1982.03320430047030> and Uno et al. (2011) <doi:10.1002/sim.4154>, the time-dependent area under the receiver operating characteristic curve (AUC) of Heagerty, Lumley and Pepe (2000) <doi:10.1111/j.0006-341X.2000.00337.x>, and the Brier score of Graf et al. (1999) <doi:10.1002/(SICI)1097-0258(19990915/30)18:17/18%3C2529::AID-SIM274%3E3.0.CO;2-5>. Supports right-censored data; competing risks, following Wolbers et al. (2014) <doi:10.1093/biostatistics/kxt059>, Wu and Li (2018) <doi:10.1002/sim.7806> and Zhuang et al. (2025) <doi:10.48550/arXiv.2507.15040>; and two-phase sampling designs such as case-cohort and nested case-control studies, following Lee, Zeleniuch-Jacquotte and Liu (2021) <doi:10.1080/02664763.2020.1861225>. Revives and expands the original 'PAmeasures' package with a unified interface for survival prediction assessment. |
| Authors: | Tong Zhu [aut], Zian Zhuang [aut, cre], Wen Su [aut], Xiaowu Dai [aut], Xiaoyan Wang [aut], Gang Li [aut] |
| Maintainer: | Zian Zhuang <[email protected]> |
| License: | GPL-3 |
| Version: | 2.0.0.9000 |
| Built: | 2026-10-29 11:58:21 UTC |
| Source: | https://github.com/larryzza/pameasures |
Transistor counts for commercial microprocessors together with their year of introduction. Moore's Law observes that the number of transistors in a dense integrated circuit roughly doubles every two years.
mooremoore
A data frame with 48 rows and 3 variables:
Year of introduction, from 1971 to 2012.
Years elapsed since 1970.
Number of transistors.
Compiled from published microprocessor transistor counts.
Assess calibration from externally supplied probabilities and observed outcomes, without requiring the original prediction model.
pa_calibration(p, event_time, status, t_star, nk = 5, case_weights = NULL) pa_calibration_cr( p, event_time, status, t_star, event_type = 1, bandwidth = NULL )pa_calibration(p, event_time, status, t_star, nk = 5, case_weights = NULL) pa_calibration_cr( p, event_time, status, t_star, event_type = 1, bandwidth = NULL )
p |
Numeric vector of predicted event probabilities at |
event_time |
Numeric vector of finite, non-negative observed times. |
status |
Event indicators: 0 for censoring and 1 for an event in
|
t_star |
One positive finite evaluation time, not exceeding the maximum
observed follow-up. Use the same time units as |
nk |
Number of restricted cubic spline knots (at least 3; default 5). |
case_weights |
Optional positive finite inverse-inclusion sampling
weights, aligned with the subjects. For |
event_type |
Positive integer code for the cause of interest. |
bandwidth |
Nearest-neighborhood smoothing bandwidth passed to
|
pa_calibration() uses pmcalibration::pmcalibration() with a
restricted cubic spline and no confidence interval computation. The knot
count is capped at the number of distinct fitting probabilities; with fewer
than five distinct probabilities or fewer than six subjects it uses the
backend's single-predictor Cox calibration instead. Constant predictions
are compared with the marginal Kaplan–Meier risk.
The right-censored backend requires probabilities strictly inside (0, 1).
Probabilities within sqrt(.Machine$double.eps) of 0 or 1 are clipped
to that distance from the boundary for fitting, with a warning (except for
the constant-prediction Kaplan–Meier branch). Metrics compare calibrated values
with the original probabilities and retain every subject.
pa_calibration_cr() uses local Aalen–Johansen estimates from
prodlim, conditional on the predicted CIF. All nonzero causes other
than event_type are pooled as competing events, not censored.
Constant or few-valued predictions use the backend's discrete risk groups.
This is a different curve estimator from Fine–Gray spline calibration.
For unweighted evaluations, both functions use
pmcalibration::cal_metrics() to summarize absolute
differences between predicted and curve-implied probabilities at the original
subjects. Only ici (the backend's Eavg), e50, and
e90 are returned. They are probability-scale errors, with 0 ideal.
ICI averages over subjects, not over follow-up time or a uniform risk grid.
With case_weights, pa_calibration() fits a sampling-weighted
Cox model directly, using the same complementary log-log transformation
and restricted cubic spline knot-placement rule as the ordinary backend.
Both coefficients and baseline cumulative hazard come from the weighted
fit. Constant predictions use weighted Kaplan–Meier risk; with fewer than
five distinct fitting probabilities or six subjects a single linear
predictor is used. Knot locations are recorded in knots.
ICI is the normalized sampling-weighted mean absolute error. E50 and E90
use weighted order-statistic interpolation, with weights normalized to the
sample size (Hmisc::wtd.quantile(..., normwt = TRUE)).
Multiplying every sampling weight by the same constant has no effect.
Supply subject inclusion weights for the actual case–cohort/NCC design, rather than weights already combined with censoring adjustments. Only sampled subjects with positive inclusion probabilities should be supplied. The calibration Cox model assumes conditionally independent censoring given its predictor and an appropriate proportional-hazards spline relationship. The competing-risks function remains an unweighted cohort evaluation. No optimism correction or design-based confidence interval is applied; preferably supply held-out or out-of-fold predictions.
A list with metrics, a named numeric vector
c(ici, e50, e90); curve, a data frame with predicted
and calibrated probabilities in input subject order;
t_star; event_type; and method, the curve estimator.
Weighted results also include the original case_weights and spline
knots on the complementary log-log scale (empty for simpler fits).
Austin PC, Steyerberg EW (2019). The Integrated Calibration Index (ICI) and related metrics for quantifying the calibration of logistic regression models. Statistics in Medicine 38, 4051–4065. doi:10.1002/sim.8281. Austin PC, Harrell FE Jr, van Klaveren D (2020). Graphical calibration curves and the integrated calibration index (ICI) for survival models. Statistics in Medicine 39, 2714–2742. doi:10.1002/sim.8570. Breslow NE, Wellner JA (2007). Weighted likelihood for semiparametric models and two-phase stratified samples, with application to Cox regression. Scandinavian Journal of Statistics 34, 86–102. doi:10.1111/j.1467-9469.2006.00523.x. Samuelsen SO (1997). A pseudolikelihood approach to analysis of nested case-control studies. Biometrika 84, 379–394. doi:10.1093/biomet/84.2.379.
set.seed(19) x <- stats::rnorm(100) event <- stats::rexp(100, rate = exp(x)) censor <- stats::rexp(100, rate = 0.3) time <- pmin(event, censor) status <- as.integer(event <= censor) p <- 1 - exp(-exp(x) * 0.5) pa_calibration(p, time, status, t_star = 0.5)$metrics # Externally supplied CIFs; no fitted model object is needed. competing <- stats::rexp(100, rate = 0.4) cr_time <- pmin(event, competing, censor) cr_status <- ifelse(censor < pmin(event, competing), 0L, ifelse(event <= competing, 1L, 2L)) cif <- exp(x) / (exp(x) + 0.4) * (1 - exp(-(exp(x) + 0.4) * 0.5)) pa_calibration_cr(cif, cr_time, cr_status, t_star = 0.5)$metricsset.seed(19) x <- stats::rnorm(100) event <- stats::rexp(100, rate = exp(x)) censor <- stats::rexp(100, rate = 0.3) time <- pmin(event, censor) status <- as.integer(event <= censor) p <- 1 - exp(-exp(x) * 0.5) pa_calibration(p, time, status, t_star = 0.5)$metrics # Externally supplied CIFs; no fitted model object is needed. competing <- stats::rexp(100, rate = 0.4) cr_time <- pmin(event, competing, censor) cr_status <- ifelse(censor < pmin(event, competing), 0L, ifelse(event <= competing, 1L, 2L)) cif <- exp(x) / (exp(x) + 0.4) * (1 - exp(-(exp(x) + 0.4) * 0.5)) pa_calibration_cr(cif, cr_time, cr_status, t_star = 0.5)$metrics
Computes static sampling weights for case-cohort evaluation. If strata
is provided, non-case weights are calibrated within each stratum.
pa_case_cohort_weights(time, status, subcohort = NULL, strata = NULL)pa_case_cohort_weights(time, status, subcohort = NULL, strata = NULL)
time |
Numeric vector of follow-up times. |
status |
Numeric {0,1}, 1 = event, 0 = censored. |
subcohort |
Logical vector; TRUE if in subcohort. |
strata |
Optional factor/character for strata (stratified cc). |
Supply the full phase-one cohort, not only sampled subjects. Cases receive weight 1. Non-cases receive the number of full-cohort non-cases divided by the number of subcohort non-cases, within each stratum if supplied. These are calibrated inverse-inclusion weights, not the time-dependent risk-set construction of the Prentice Cox estimator. Passing these as ordinary Cox case weights does not implement that estimator.
Retain only the random subcohort and all cases for evaluation, and align weights with the evaluation subjects. At least one sampled non-case is needed in every stratum containing non-cases to obtain finite weights.
Numeric vector of weights in full-cohort input order. Weights for non-sampled non-cases are also returned and must be excluded from the sampled analysis dataset.
# Compute weights using the full cohort, then retain the random subcohort # and all cases for analysis. Sampled non-cases represent cohort non-cases. d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) set.seed(2025) w <- pa_case_cohort_weights(time = d$time, status = d$status, subcohort = stats::rbinom(nrow(d), 1, 0.4)) summary(w)# Compute weights using the full cohort, then retain the random subcohort # and all cases for analysis. Sampled non-cases represent cohort non-cases. d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) set.seed(2025) w <- pa_case_cohort_weights(time = d$time, status = d$status, subcohort = stats::rbinom(nrow(d), 1, 0.4)) summary(w)
This function constructs a standardized evaluation dataset (dat1) from the
analysis dataset, predicted survival probabilities, and sampling weights, and
computes key performance metrics in a tidy table format.
Works for common two-phase sampling designs:
- Case-cohort (unstratified or stratified)
- Nested case-control (matched or unmatched)
pa_evaluate_two_phase( pred_results, t_star = NULL, tau = 1e+11, km_cens_fit, case_weights, metrics = c("pseudo_r2", "harrell_c", "uno_c", "brier_score", "td_auc", "ici"), calibration_args = list() )pa_evaluate_two_phase( pred_results, t_star = NULL, tau = 1e+11, km_cens_fit, case_weights, metrics = c("pseudo_r2", "harrell_c", "uno_c", "brier_score", "td_auc", "ici"), calibration_args = list() )
pred_results |
A prediction list, such as the output of
|
t_star |
(Optional) Evaluation time for Brier score, time-dependent AUC, and calibration. If NULL, it defaults to the median observed follow-up time. |
tau |
(Optional) A numeric value specifying the truncation time for calculating explained variation metrics. If NULL, it defaults to the maximum observed survival time. |
km_cens_fit |
Kaplan-Meier estimate of the censoring distribution, a
|
case_weights |
Numeric sampling weights in the same subject order as
|
metrics |
Character vector of metrics. |
calibration_args |
Optional named list of calibration settings,
e.g. |
Calibration uses the event risk from the last prediction-grid
column at or before t_star, with the exact requested outcome
horizon. Both the Cox spline calibration model and absolute-error
summaries use sampling weights; see pa_calibration().
Calibration-only inputs need times, status, and
surv_prob, without fitted models or mean time predictions.
A tibble with columns:
- Metric: name of metric
- Value: numeric value
d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # case-cohort sampling weights, and the censoring distribution they need set.seed(2025) in_subcohort <- stats::rbinom(nrow(d), 1, 0.4) == 1 w <- pa_case_cohort_weights(d$time, d$status, subcohort = in_subcohort) keep <- which(in_subcohort | d$status == 1) keep <- keep[order(d$time[keep])] pred <- pa_predict_coxph(fit, c("x1", "x2"), new_data = d[keep, ]) w <- w[keep] km <- survival::survfit(survival::Surv(d$time, 1 - d$status) ~ 1) pa_evaluate_two_phase(pred_results = pred, km_cens_fit = km, case_weights = w)d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # case-cohort sampling weights, and the censoring distribution they need set.seed(2025) in_subcohort <- stats::rbinom(nrow(d), 1, 0.4) == 1 w <- pa_case_cohort_weights(d$time, d$status, subcohort = in_subcohort) keep <- which(in_subcohort | d$status == 1) keep <- keep[order(d$time[keep])] pred <- pa_predict_coxph(fit, c("x1", "x2"), new_data = d[keep, ]) w <- w[keep] km <- survival::survfit(survival::Surv(d$time, 1 - d$status) ~ 1) pa_evaluate_two_phase(pred_results = pred, km_cens_fit = km, case_weights = w)
This function computes a comprehensive set of performance metrics for survival analysis models. It provides metrics such as R_square, L_square, Pseudo_R, Harrell's C, Uno's C, Brier Score, and Time-dependent AUC. Users can specify particular metrics and model types, enabling tailored performance evaluation for various survival models.
pa_fit_and_eval( train_data, covariates, models = "coxph", metrics = NULL, predicted_data = NULL, t_star = NULL, tau = NULL, calibration_args = list() )pa_fit_and_eval( train_data, covariates, models = "coxph", metrics = NULL, predicted_data = NULL, t_star = NULL, tau = NULL, calibration_args = list() )
train_data |
A data frame containing the survival data. |
covariates |
A character vector of covariate names to include in the model. |
models |
Character vector of model types: "coxph", "exp", "lognormal", or "weibull". The default fits only "coxph"; use "all" to fit all four. |
metrics |
A character string or vector specifying the metrics to compute.
The default
|
predicted_data |
(Optional) A data frame containing validation data. If |
t_star |
(Optional) Positive evaluation time for Brier score, AUC, and calibration; the median evaluation follow-up if omitted. |
tau |
(Optional) A time point for truncating the survival time. If provided, the function evaluates predictions up to this time point. |
calibration_args |
Optional named list passed to |
Calibration evaluates model-predicted event probabilities exactly
at t_star (the median evaluation follow-up if omitted).
Rows with missing values in time, status, or a
requested covariate are excluded from training and evaluation with a
warning. All requested metrics use the same complete evaluation rows.
The warning can be silenced with suppressWarnings().
A data frame containing the selected model's performance metrics.
data(pbc, package = "survival") pbc <- pbc[!is.na(pbc$trt), ] pbc$log_albumin <- log(pbc$albumin) pbc$log_bili <- log(pbc$bili) pbc$log_protime <- log(pbc$protime) pbc$status <- ifelse(pbc$status == 2, 1, 0) covariates <- c("age", "log_albumin", "log_bili", "log_protime", "edema") dat <- pbc[, c("time", "status", covariates)] dat <- dat[stats::complete.cases(dat), ] # Default Cox model and metrics, including ICI results <- pa_fit_and_eval(train_data = dat, covariates = covariates) results # A specific subset of models and metrics results2 <- pa_fit_and_eval( train_data = dat, covariates = covariates, models = c("lognormal", "weibull"), metrics = c("r_square", "l_square", "brier_score") ) results2data(pbc, package = "survival") pbc <- pbc[!is.na(pbc$trt), ] pbc$log_albumin <- log(pbc$albumin) pbc$log_bili <- log(pbc$bili) pbc$log_protime <- log(pbc$protime) pbc$status <- ifelse(pbc$status == 2, 1, 0) covariates <- c("age", "log_albumin", "log_bili", "log_protime", "edema") dat <- pbc[, c("time", "status", covariates)] dat <- dat[stats::complete.cases(dat), ] # Default Cox model and metrics, including ICI results <- pa_fit_and_eval(train_data = dat, covariates = covariates) results # A specific subset of models and metrics results2 <- pa_fit_and_eval( train_data = dat, covariates = covariates, models = c("lognormal", "weibull"), metrics = c("r_square", "l_square", "brier_score") ) results2
The metric identifiers accepted by the metrics argument of the evaluation
functions, and emitted in the Metric column of their results.
pa_metric_names()pa_metric_names()
A character vector of canonical metric names.
pa_metric_names()pa_metric_names()
Computes inverse-selection weights for NCC designs. If strata is
provided, matched NCC weights are computed.
pa_nested_case_control_weights(time, status, strata = NULL, m = NULL)pa_nested_case_control_weights(time, status, strata = NULL, m = NULL)
time |
Numeric vector of follow-up times. |
status |
Numeric {0,1}, 1 = event, 0 = censored. |
strata |
Optional factor/character defining matching strata in the full phase-one cohort, not IDs of sampled case-control sets. |
m |
Positive integer, number of controls per case (required). |
Supply full-cohort follow-up and event information, not only sampled rows.
Controls are sampled without replacement within each case's risk set and
independently across cases. At tied failure times, all simultaneous cases
are excluded from the eligible controls. Events precede censoring at ties.
Matching restricts both cases and controls to the same stratum. If fewer
than m controls are eligible, all eligible controls are sampled.
A non-case with no sampling opportunity has infinite inverse probability;
such a subject cannot belong to the sampled analysis dataset. Subset the
returned weights to the subjects actually sampled and align their order
with the prediction objects before evaluation.
Numeric vector of weights, in full-cohort input order.
# Unmatched nested case-control design with m = 2 controls per case. d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) set.seed(2025) w <- pa_nested_case_control_weights(time = d$time, status = d$status, m = 2) summary(w)# Unmatched nested case-control design with m = 2 controls per case. d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) set.seed(2025) w <- pa_nested_case_control_weights(time = d$time, status = d$status, m = 2) summary(w)
Evaluate observed and predicted continuous responses using the R-squared and L-squared measures of Li and Wang (2019). Predictions may come from a linear or nonlinear model; the original fitted model is not required.
pa_nlm(y, y.predict) pam.nlm(y, y.predict)pa_nlm(y, y.predict) pam.nlm(y, y.predict)
y |
A finite numeric vector of observed continuous responses. |
y.predict |
A finite numeric vector of predictions, paired with
|
An ordinary least-squares regression of y on
y.predict, with an intercept, gives the corrected predictions
. The measures are
and
R-squared describes explained response variation after affine correction;
L-squared compares the corrected and original prediction errors. These
are the original PAmeasures definitions, not .
Constant predictions use the intercept-only correction mean(y).
A constant observed response makes R-squared undefined. Predictions
exactly equal to the observations make L-squared undefined because its
denominator is zero. Undefined measures are NA_real_, with a
warning; any other defined measure is retained. Inputs whose squared
deviations overflow or whose nonzero deviations underflow to zero are
rejected with an error and should be rescaled to common units.
pam.nlm() is the retained PAmeasures 0.1.0 compatibility entry
point. It uses the same calculation, then rounds and formats each result
to four decimal places as in the historical package. Its ordinary,
defined results retain the historical field names and character format;
input validation and undefined-result handling follow pa_nlm().
No deprecation warning is issued for this supported compatibility entry.
A list with R.squared and L.squared.
pa_nlm() returns unrounded numeric scalars.
pam.nlm() returns the historical formatted character scalars.
Li, G. and Wang, X. (2019). Prediction Accuracy Measures for a Nonlinear Model and for Right-Censored Time-to-Event Data. Journal of the American Statistical Association, 114, 1815–1825. doi:10.1080/01621459.2018.1515079.
y <- c(1, 2, 4, 5) prediction <- c(1, 2, 3, 4) pa_nlm(y, prediction) pam.nlm(y, prediction) data(moore) fit <- stats::lm(log2(count) ~ time, data = moore) pa_nlm(moore$count, 2^stats::predict(fit))y <- c(1, 2, 4, 5) prediction <- c(1, 2, 3, 4) pa_nlm(y, prediction) pam.nlm(y, prediction) data(moore) fit <- stats::lm(log2(count) ~ time, data = moore) pa_nlm(moore$count, 2^stats::predict(fit))
Plots a model-derived survival-time summary and observed follow-up times against the model's linear predictor. Event and censoring status are shown with different point shapes. This is a descriptive display, not a probability-calibration curve.
pa_plot_pred( data, title = NULL, xlab = "Risk Score", ylab = "Days", label_name = c("censored", "event"), levels = NULL, shape_style = NULL, legend_name = "status", invert_linear = TRUE, sample_index = NULL, restrict_time = NULL )pa_plot_pred( data, title = NULL, xlab = "Risk Score", ylab = "Days", label_name = c("censored", "event"), levels = NULL, shape_style = NULL, legend_name = "status", invert_linear = TRUE, sample_index = NULL, restrict_time = NULL )
data |
A data frame (or list-like object) containing at least the following columns:
|
title |
Optional plot title (default = |
xlab, ylab
|
Character strings specifying the axis labels
(defaults are |
label_name |
Character vector of labels corresponding to event
status codes (e.g., |
levels |
Numeric or character vector giving the ordering of status
levels; defaults inferred from the length of |
shape_style |
Numeric vector of plotting symbols (pch values)
corresponding to the levels of |
legend_name |
Character string for the legend title (default = "status"). |
invert_linear |
Logical. If |
sample_index |
Optional integer vector selecting which subjects to
plot, e.g. |
restrict_time |
Optional numeric. If supplied, observed times above this value are capped at it and a dashed horizontal reference line is drawn, which keeps a few long follow-up times from compressing the plot. |
The function is primarily designed for descriptive diagnostics in survival
model development. The line (pred vs linear.pred) shows the
predicted values joined in x order, without smoothing, while the overlaid points (times vs
linear.pred) illustrate observed event times, distinguished by
censoring status. For more elaborate visualization (e.g., confidence bands,
multiple models), use faceting or combine plots with patchwork.
A ggplot object that can be further customized.
d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # observed time against predicted risk score, shaped by event status pa_plot_pred(pred, title = "Cox model", invert_linear = FALSE)d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # observed time against predicted risk score, shaped by event status pa_plot_pred(pred, title = "Cox model", invert_linear = FALSE)
Builds one panel per data frame using plot_fun (default: pa_plot_pred)
and arranges them with letter tags. The function collects a single shared legend
and applies global styling (titles, axis labels, legend options) across panels.
pa_plot_summary( data_list, titles = NULL, plot_fun = pa_plot_pred, ncol = 2, tag_levels = "a", invert_linear = TRUE, xlab = "Risk Score", ylab = "Days", label_name = c("Censored", "Event"), levels = NULL, shape_style = NULL, sample_size = NULL, restrict_time = NULL, legend_name = "Status", seed = NULL )pa_plot_summary( data_list, titles = NULL, plot_fun = pa_plot_pred, ncol = 2, tag_levels = "a", invert_linear = TRUE, xlab = "Risk Score", ylab = "Days", label_name = c("Censored", "Event"), levels = NULL, shape_style = NULL, sample_size = NULL, restrict_time = NULL, legend_name = "Status", seed = NULL )
data_list |
List of data frames, each containing columns:
|
titles |
Optional character vector, same length as |
plot_fun |
Function that creates a single panel (default = |
ncol |
Number of columns in the panel layout (default = 2). |
tag_levels |
Letter style for tags: |
invert_linear |
Logical or logical vector; if length 1, recycled to all
panels. Use |
xlab, ylab
|
Axis labels applied to all panels (defaults are "Risk Score" and "Days"). |
label_name |
Character vector of labels for status legend (e.g., |
levels |
Numeric vector for |
shape_style |
Numeric vector of plotting symbols for censoring/event status. |
sample_size |
Optional integer. If supplied, that many subjects are
sampled at random for each panel using |
restrict_time |
Optional numeric. If supplied, observed times above this value are capped at it and a dashed horizontal reference line is drawn, which keeps a few long follow-up times from compressing the plot. |
legend_name |
Character string for the legend title (default = "status"). |
seed |
Optional integer. If supplied, the random number generator is
seeded with it before |
A patchwork ggplot object combining all panels.
# One panel per model, with larger x values indicating higher risk. d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] cox <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) wei <- survival::survreg(survival::Surv(time, status) ~ x1 + x2, data = d, dist = "weibull", x = TRUE, y = TRUE) p_cox <- pa_predict_coxph(model = cox, covs = c("x1", "x2"), new_data = d) p_wei <- pa_predict_survreg(model = wei, covs = c("x1", "x2"), new_data = d) pa_plot_summary( list(p_cox, p_wei), titles = c("Cox PH", "Weibull AFT"), invert_linear = c(FALSE, TRUE), ncol = 2, legend_name = "Status" )# One panel per model, with larger x values indicating higher risk. d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] cox <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) wei <- survival::survreg(survival::Surv(time, status) ~ x1 + x2, data = d, dist = "weibull", x = TRUE, y = TRUE) p_cox <- pa_predict_coxph(model = cox, covs = c("x1", "x2"), new_data = d) p_wei <- pa_predict_survreg(model = wei, covs = c("x1", "x2"), new_data = d) pa_plot_summary( list(p_cox, p_wei), titles = c("Cox PH", "Weibull AFT"), invert_linear = c(FALSE, TRUE), ncol = 2, legend_name = "Status" )
This function calculates cumulative incidence function (CIF) predictions for individual subjects across multiple types of competing risks models. It supports cause-specific Cox models, accelerated failure time (AFT) models (Weibull, Log-normal, Log-logistic), Fine-Gray models, and competing risks random survival forests.
pa_predict_cif( model1 = NULL, model2 = NULL, fg_model = NULL, cr_model = NULL, tau = NULL, newdata, event.type = 1, covs )pa_predict_cif( model1 = NULL, model2 = NULL, fg_model = NULL, cr_model = NULL, tau = NULL, newdata, event.type = 1, covs )
model1 |
The fitted model for the target cause (either a |
model2 |
The fitted model for competing causes (either a |
fg_model |
Optional Fine-Gray model object ( |
cr_model |
Optional competing risks random forest model object. |
tau |
The maximum time horizon for prediction. |
newdata |
A data frame containing observed event times ( |
event.type |
The event type of interest for which CIF is to be predicted. |
covs |
A character vector specifying the names of the covariates used in the model. |
For cause-specific AFT models, CIFs are evaluated by adaptive quadrature on the target cause's probability scale. This handles Weibull densities that are unbounded at time zero, including fixed scales above 1. The time summary integrates one minus the CIF on the returned time grid.
A list with the following components:
Numeric vector of observed times from newdata$time.
Integer vector of event indicators from newdata$status.
A matrix where the first column is time, and subsequent columns are predicted CIFs for each subject.
Numeric vector of the integral of one minus the CIF from 0 to tau.
Numeric vector of subject-specific linear predictors.
data(pbc, package = "survival") pbc <- pbc[!is.na(pbc$trt), ] pbc$log_albumin <- log(pbc$albumin) pbc$log_bili <- log(pbc$bili) pbc$log_protime <- log(pbc$protime) covariates <- c("age", "log_albumin", "log_bili", "log_protime", "edema") dat <- pbc[, c("time", "status", covariates)] dat <- dat[stats::complete.cases(dat), ] # Cause-specific Cox models: event of interest (1) and competing event (2) m1 <- survival::coxph(survival::Surv(time, status == 1) ~ age + log_albumin + log_bili + log_protime + edema, data = dat, x = TRUE, y = TRUE) m2 <- survival::coxph(survival::Surv(time, status == 2) ~ age + log_albumin + log_bili + log_protime + edema, data = dat, x = TRUE, y = TRUE) pred <- pa_predict_cif(model1 = m1, model2 = m2, newdata = dat, covs = covariates, event.type = 1, tau = max(dat$time)) str(pred$pred) # Fine-Gray backend, available when cmprsk is installed if (requireNamespace("cmprsk", quietly = TRUE)) { fg <- cmprsk::crr(ftime = dat$time, fstatus = dat$status, cov1 = as.matrix(dat[, covariates]), failcode = 1, cencode = 0) pred_fg <- pa_predict_cif(fg_model = fg, newdata = dat, covs = covariates, event.type = 1, tau = max(dat$time)) str(pred_fg$pred) }data(pbc, package = "survival") pbc <- pbc[!is.na(pbc$trt), ] pbc$log_albumin <- log(pbc$albumin) pbc$log_bili <- log(pbc$bili) pbc$log_protime <- log(pbc$protime) covariates <- c("age", "log_albumin", "log_bili", "log_protime", "edema") dat <- pbc[, c("time", "status", covariates)] dat <- dat[stats::complete.cases(dat), ] # Cause-specific Cox models: event of interest (1) and competing event (2) m1 <- survival::coxph(survival::Surv(time, status == 1) ~ age + log_albumin + log_bili + log_protime + edema, data = dat, x = TRUE, y = TRUE) m2 <- survival::coxph(survival::Surv(time, status == 2) ~ age + log_albumin + log_bili + log_protime + edema, data = dat, x = TRUE, y = TRUE) pred <- pa_predict_cif(model1 = m1, model2 = m2, newdata = dat, covs = covariates, event.type = 1, tau = max(dat$time)) str(pred$pred) # Fine-Gray backend, available when cmprsk is installed if (requireNamespace("cmprsk", quietly = TRUE)) { fg <- cmprsk::crr(ftime = dat$time, fstatus = dat$status, cov1 = as.matrix(dat[, covariates]), failcode = 1, cencode = 0) pred_fg <- pa_predict_cif(fg_model = fg, newdata = dat, covs = covariates, event.type = 1, tau = max(dat$time)) str(pred_fg$pred) }
Returns subject-specific survival probabilities and a time summary from a
fitted model when predict = TRUE. With predict = FALSE, returns
R-squared, L-squared and pseudo-R-squared instead. R-squared measures explained
variation after affine correction. L-squared is the corrected prediction
error divided by the original prediction error; it is not the fraction of
error removed by correction. Pseudo-R-squared is their product.
pa_predict_coxph(model, covs, tau = 1e+11, new_data = NULL, predict = TRUE)pa_predict_coxph(model, covs, tau = 1e+11, new_data = NULL, predict = TRUE)
model |
An object of class |
covs |
A character vector specifying the names of the covariates used in the model. |
tau |
(Optional) A restriction time for computing restricted survival time. If provided, the function evaluates predictions up to this time point. Default is 10e10. |
new_data |
Optional evaluation data. Its outcomes are returned for metric computation and define the evaluation time grid; they never estimate the baseline hazard. All prediction parameters come from the training fit. |
predict |
A logical value indicating whether to return individual predictions. Default is TRUE. |
Probabilities use the fitted training baseline, including its case
weights, strata, tie method and offsets. Offset centering follows the
weighted training mean used by survfit().
With predict = FALSE, R-squared metrics return missing values with a
warning if censoring weights cannot be normalized, for example when the
training censoring survival is zero at an observed evaluation event.
With predict = TRUE, a single prediction list containing:
model – the fitted survival model object.
times – numeric vector of observed follow-up times for the evaluation dataset.
status – event indicators (1 = event, 0 = censored).
surv_prob – an matrix of predicted
subject-specific survival probabilities on a common time grid.
pred – numerical integral of the training-model survival curve
on the evaluation grid, up to its maximum time or tau, whichever is
smaller. This preserves the legacy integration rule; it is not an
identified unrestricted mean. Fix the grid and tau for comparisons.
linear.pred – the model's linear predictor, aligned with
the returned subjects.
covs – character vector of covariate names used for prediction.
new_data – dataset on which the predictions were computed.
Subjects are sorted by observed follow-up time. The columns of
surv_prob correspond to the sorted follow-up times, including ties.
With predict = FALSE, the result is a list of formatted character
scalars named R.squared, L.squared, and Psuedo.R.
The last spelling is retained for compatibility.
Li, G., & Wang, X. (2016). Prediction Accuracy Measures for a Nonlinear Model and for Right-Censored Time-to-Event Data. arXiv preprint arXiv:1611.03063. Available at https://arxiv.org/abs/1611.03063.
data(pbc, package = "survival") pbc <- pbc[!is.na(pbc$trt), ] pbc$log_albumin <- log(pbc$albumin) pbc$log_bili <- log(pbc$bili) pbc$log_protime <- log(pbc$protime) pbc$status <- ifelse(pbc$status == 2, 1, 0) covs <- c("age", "log_albumin", "log_bili", "log_protime", "edema") df <- pbc[, c("time", "status", covs)] df <- df[stats::complete.cases(df), ] set.seed(2025) idx <- sample(seq_len(nrow(df)), round(2 / 3 * nrow(df))) train_data <- df[idx, ] test_data <- df[-idx, ] m.cox <- survival::coxph(survival::Surv(time, status) ~ ., data = train_data, x = TRUE, y = TRUE) cox_pred <- pa_predict_coxph(model = m.cox, covs = covs, new_data = test_data) str(cox_pred$pred)data(pbc, package = "survival") pbc <- pbc[!is.na(pbc$trt), ] pbc$log_albumin <- log(pbc$albumin) pbc$log_bili <- log(pbc$bili) pbc$log_protime <- log(pbc$protime) pbc$status <- ifelse(pbc$status == 2, 1, 0) covs <- c("age", "log_albumin", "log_bili", "log_protime", "edema") df <- pbc[, c("time", "status", covs)] df <- df[stats::complete.cases(df), ] set.seed(2025) idx <- sample(seq_len(nrow(df)), round(2 / 3 * nrow(df))) train_data <- df[idx, ] test_data <- df[-idx, ] m.cox <- survival::coxph(survival::Surv(time, status) ~ ., data = train_data, x = TRUE, y = TRUE) cox_pred <- pa_predict_coxph(model = m.cox, covs = covs, new_data = test_data) str(cox_pred$pred)
Returns subject-specific survival probabilities and a time summary from a
fitted model when predict = TRUE. With predict = FALSE, returns
R-squared, L-squared and pseudo-R-squared instead. R-squared measures explained
variation after affine correction. L-squared is the corrected prediction
error divided by the original prediction error; it is not the fraction of
error removed by correction. Pseudo-R-squared is their product.
pa_predict_survreg(model, covs, tau = 1e+11, new_data = NULL, predict = TRUE)pa_predict_survreg(model, covs, tau = 1e+11, new_data = NULL, predict = TRUE)
model |
An object of class |
covs |
A vector of covariate names to be used for prediction when new_data is provided. |
tau |
(Optional) A restriction time for computing restricted survival time. If provided, the function evaluates predictions up to this time point. Default is 10e10. |
new_data |
(Optional) A new dataset for evaluating the model. If NULL, the function uses the training data stored in model. |
predict |
(Optional) A logical value indicating whether to return individual predictions. Default is TRUE. |
With predict = FALSE, R-squared metrics return missing values
with a warning if censoring weights cannot be normalized, for example
when the training censoring survival is zero at an observed evaluation
event. Individual predictions remain available with predict = TRUE.
With predict = TRUE, a single prediction list containing:
model – the fitted survival model object.
times – numeric vector of observed follow-up times for the evaluation dataset.
status – event indicators (1 = event, 0 = censored).
surv_prob – an matrix of predicted
subject-specific survival probabilities on a common time grid.
pred – model-based restricted mean survival time
, calculated analytically up to
tau; it is not the median survival time.
linear.pred – the model's linear predictor, aligned with
the returned subjects.
covs – character vector of covariate names used for prediction.
new_data – dataset on which the predictions were computed.
Subjects are sorted by observed follow-up time. The columns of
surv_prob correspond to the sorted follow-up times, including ties.
With predict = FALSE, the result is a list of formatted character
scalars named R.squared, L.squared, and Psuedo.R.
The last spelling is retained for compatibility.
data(pbc, package = "survival") pbc <- pbc[!is.na(pbc$trt), ] pbc$log_albumin <- log(pbc$albumin) pbc$log_bili <- log(pbc$bili) pbc$log_protime <- log(pbc$protime) pbc$status <- ifelse(pbc$status == 2, 1, 0) covs <- c("age", "log_albumin", "log_bili", "log_protime", "edema") df <- pbc[, c("time", "status", covs)] df <- df[stats::complete.cases(df), ] set.seed(2025) idx <- sample(seq_len(nrow(df)), round(2 / 3 * nrow(df))) train_data <- df[idx, ] test_data <- df[-idx, ] m.wei <- survival::survreg(survival::Surv(time, status) ~ ., data = train_data, dist = "weibull", x = TRUE, y = TRUE) wei_pred <- pa_predict_survreg(model = m.wei, covs = covs, new_data = test_data) str(wei_pred$pred)data(pbc, package = "survival") pbc <- pbc[!is.na(pbc$trt), ] pbc$log_albumin <- log(pbc$albumin) pbc$log_bili <- log(pbc$bili) pbc$log_protime <- log(pbc$protime) pbc$status <- ifelse(pbc$status == 2, 1, 0) covs <- c("age", "log_albumin", "log_bili", "log_protime", "edema") df <- pbc[, c("time", "status", covs)] df <- df[stats::complete.cases(df), ] set.seed(2025) idx <- sample(seq_len(nrow(df)), round(2 / 3 * nrow(df))) train_data <- df[idx, ] test_data <- df[-idx, ] m.wei <- survival::survreg(survival::Surv(time, status) ~ ., data = train_data, dist = "weibull", x = TRUE, y = TRUE) wei_pred <- pa_predict_survreg(model = m.wei, covs = covs, new_data = test_data) str(wei_pred$pred)
For a list of models (each providing times, status, surv_prob, pred),
calls pa_evaluate_two_phase() for each and returns a wide
comparison table with metrics as rows and model names as columns.
pa_sample_design( models, case_weights, km_cens, metrics = c("pseudo_r2", "harrell_c", "uno_c", "brier_score", "td_auc", "ici"), t_star = NULL, tau = NULL, digits = 2, calibration_args = list() )pa_sample_design( models, case_weights, km_cens, metrics = c("pseudo_r2", "harrell_c", "uno_c", "brier_score", "td_auc", "ici"), t_star = NULL, tau = NULL, digits = 2, calibration_args = list() )
models |
A named list of prediction objects, such as those returned by
|
case_weights |
Numeric vector of inverse-inclusion sampling weights, aligned with each model's evaluation subjects. |
km_cens |
A Kaplan-Meier fit for censoring (based on training data) used to compute IPCW. Not required for calibration-only evaluation. |
metrics |
Character vector of metrics to compute (default shown below). |
t_star |
Optional evaluation time for Brier score, AUC and calibration. Defaults to the median observed follow-up time. |
tau |
Optional truncation time for explained-variation metrics.
If |
digits |
Integer number of decimal places to round values (default 2). |
calibration_args |
Optional named list passed to |
A data.frame: rows = metrics, columns = model names.
d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # case-cohort sampling weights, and the censoring distribution they need set.seed(2025) in_subcohort <- stats::rbinom(nrow(d), 1, 0.4) == 1 w <- pa_case_cohort_weights(d$time, d$status, subcohort = in_subcohort) keep <- which(in_subcohort | d$status == 1) keep <- keep[order(d$time[keep])] pred <- pa_predict_coxph(fit, c("x1", "x2"), new_data = d[keep, ]) w <- w[keep] km <- survival::survfit(survival::Surv(d$time, 1 - d$status) ~ 1) # the wide summary of the same two-phase evaluation pa_sample_design(models = list(cc = pred), case_weights = w, km_cens = km)d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # case-cohort sampling weights, and the censoring distribution they need set.seed(2025) in_subcohort <- stats::rbinom(nrow(d), 1, 0.4) == 1 w <- pa_case_cohort_weights(d$time, d$status, subcohort = in_subcohort) keep <- which(in_subcohort | d$status == 1) keep <- keep[order(d$time[keep])] pred <- pa_predict_coxph(fit, c("x1", "x2"), new_data = d[keep, ]) w <- w[keep] km <- survival::survfit(survival::Surv(d$time, 1 - d$status) ~ 1) # the wide summary of the same two-phase evaluation pa_sample_design(models = list(cc = pred), case_weights = w, km_cens = km)
Generates survival data under a Cox proportional hazards model with a Weibull
baseline cumulative hazard . Event times are sampled via
inverse transform:
where .
Covariates are simulated from a standard normal distribution.
With both interact = FALSE and nonlinear = FALSE, the design
matrix contains only main effects.
When interact = TRUE, the function automatically detects whether length(beta)
equals 3 or 7:
If length(beta) = 3, two covariates () and their
two-way interaction () are included.
If length(beta) = 7, three covariates (),
their three two-way interactions (),
and their three-way interaction () are included.
Right-censoring times are generated as with
. By default, sd is set to
, and mu is calibrated so that the expected censoring
proportion equals pi_c. If mu and/or sd are supplied, calibration
is skipped accordingly:
mu and sd provided: use them directly (no calibration).
sd only: solve for mu using the given sd.
mu only: use the given mu and estimate .
neither: estimate both and mu.
pa_sim_cox_weibull( n, pi_c, v, beta, mu = NULL, sd = NULL, seed = NULL, interact = FALSE, nonlinear = FALSE )pa_sim_cox_weibull( n, pi_c, v, beta, mu = NULL, sd = NULL, seed = NULL, interact = FALSE, nonlinear = FALSE )
n |
Integer. Sample size. |
pi_c |
Numeric in |
v |
Positive numeric. Weibull shape in |
beta |
Numeric vector of regression coefficients.
If |
mu |
Optional numeric. Mean of |
sd |
Optional positive numeric. Standard deviation used in
|
seed |
Optional integer seed for reproducibility. |
interact |
Logical (default |
nonlinear |
Logical. If |
A data.frame with columns:
time: observed time .
status: event indicator .
x1, x2, ...: simulated covariates.
y_true: latent event time , only when pi_c > 0.
cens_time: censoring time , only when pi_c > 0.
Interaction and quadratic terms are used internally to generate event
times, but only the original covariates are returned as columns.
With pi_c = 0, time is the latent event time and all
status values are 1; no model or censoring attributes are attached.
With pi_c > 0, additional attributes are:
mu: censoring shift used.
sd_log: standard deviation used in .
sd_logY: from the generated events.
pi_c_target: requested censoring rate.
pi_c_observed: realized censoring rate.
v, beta: model parameters.
# (1) Main effects only set.seed(1) d0 <- pa_sim_cox_weibull(500, 0.3, v = 1.2, beta = c(0.4, -0.2)) # (2) Two covariates + their interaction d1 <- pa_sim_cox_weibull(500, 0.3, v = 1.2, beta = c(0.5, -0.3, 0.2), interact = TRUE) # (3) Three covariates + all two-way and three-way interactions b <- c(0.4, -0.2, 0.1, 0.15, -0.05, 0.08, 0.02) d2 <- pa_sim_cox_weibull(500, 0.3, v = 1.0, beta = b, interact = TRUE)# (1) Main effects only set.seed(1) d0 <- pa_sim_cox_weibull(500, 0.3, v = 1.2, beta = c(0.4, -0.2)) # (2) Two covariates + their interaction d1 <- pa_sim_cox_weibull(500, 0.3, v = 1.2, beta = c(0.5, -0.3, 0.2), interact = TRUE) # (3) Three covariates + all two-way and three-way interactions b <- c(0.4, -0.2, 0.1, 0.15, -0.05, 0.08, 0.02) d2 <- pa_sim_cox_weibull(500, 0.3, v = 1.0, beta = b, interact = TRUE)
This function generates a simulated dataset for a competing risks model with two causes, following a Fine-Gray-like model specification.
pa_simulate_fine_gray( n, v, beta1, beta2, lambda1 = 1, X = NULL, mu = 0, p = 0.7, c_scale = 1, censor = 0, sd.time = NULL, mu.c.time = NULL, independent_c = TRUE, report.mu_and_sd = FALSE, seed = NULL )pa_simulate_fine_gray( n, v, beta1, beta2, lambda1 = 1, X = NULL, mu = 0, p = 0.7, c_scale = 1, censor = 0, sd.time = NULL, mu.c.time = NULL, independent_c = TRUE, report.mu_and_sd = FALSE, seed = NULL )
n |
An integer specifying the number of subjects to simulate. |
v |
A numeric value for the Weibull shape parameter. |
beta1 |
A numeric vector of coefficients for cause 1. |
beta2 |
A numeric vector of coefficients for cause 2. |
lambda1 |
Positive coefficient in the baseline term
|
X |
A matrix of covariates, if not provided, covariates are generated from a standard normal distribution. |
mu |
A numeric value for the mean of the generated covariates. |
p |
Baseline probability of a cause 1 event when the cause-1 linear
predictor is zero. Conditional on covariates, this probability is
|
c_scale |
Positive multiplier for the default standard deviation of log censoring times. |
censor |
Target censoring proportion. Zero produces uncensored data.
For a positive value, the log-censoring shift is calibrated to this target
unless |
sd.time |
Standard deviation of log censoring times. If |
mu.c.time |
Optional additive shift of log censoring times. When
|
independent_c |
A logical value. If TRUE, censoring times are independent of covariates. |
report.mu_and_sd |
A logical value. If TRUE, a summary of censoring parameters and event time quantiles is returned instead of the full dataset. |
seed |
Optional integer. If supplied, the random number generator is
seeded with it for reproducibility. The default |
A data frame containing the simulated data with columns obs.times, obs.event, and the covariates X.
Alternatively, if report.mu_and_sd is TRUE, it returns a data frame with censoring
parameters and event time quantiles.
# Simulate a dataset with 100 subjects sim_data <- pa_simulate_fine_gray( n = 100, v = 1, beta1 = c(0.5, -0.2), beta2 = c(-0.3, 0.4) ) head(sim_data) # Simulate with a specific censoring proportion and report parameters params_summary <- pa_simulate_fine_gray( n = 1000, v = 1, beta1 = c(0.5, -0.2), beta2 = c(-0.3, 0.4), censor = 0.2, report.mu_and_sd = TRUE ) params_summary# Simulate a dataset with 100 subjects sim_data <- pa_simulate_fine_gray( n = 100, v = 1, beta1 = c(0.5, -0.2), beta2 = c(-0.3, 0.4) ) head(sim_data) # Simulate with a specific censoring proportion and report parameters params_summary <- pa_simulate_fine_gray( n = 1000, v = 1, beta1 = c(0.5, -0.2), beta2 = c(-0.3, 0.4), censor = 0.2, report.mu_and_sd = TRUE ) params_summary
Calls pa_survival_eval() for each model in a named list and
outputs a wide table: each row corresponds to one metric and each column
corresponds to a model.
pa_summarize( models, metrics = NULL, t_star = NULL, tau = 1e+11, digits = 2, calibration_args = list() )pa_summarize( models, metrics = NULL, t_star = NULL, tau = 1e+11, digits = 2, calibration_args = list() )
models |
A named list in which each element corresponds to a fitted model. Each model entry must itself be a list containing:
Note: If |
metrics |
Optional character vector of metrics to compute (passed through). |
t_star |
Optional evaluation time for point metrics and calibration. Defaults to the median observed follow-up. |
tau |
Optional numeric scalar, specify the max time horizon for R2 measure and concordance measure (default = 10e10). |
digits |
Number of decimal places to round the metric values (default = 2). |
calibration_args |
Optional named list passed to |
A data frame with:
Each row = metric
Each column = model name
Cell values = corresponding rounded metric values
d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # one column per model, metrics down the rows pa_summarize(list(cox = pred))d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # one column per model, metrics down the rows pa_summarize(list(cox = pred))
Calls pa_survival_eval_cr() for each model in a named list and
returns a wide table where rows are metrics and columns are model names.
pa_summarize_cr( models, metrics = NULL, t_star = NULL, tau = NULL, event_type = 1, digits = 2, censoring = c("cox", "km"), censor_covariates = NULL, calibration_args = list() )pa_summarize_cr( models, metrics = NULL, t_star = NULL, tau = NULL, event_type = 1, digits = 2, censoring = c("cox", "km"), censor_covariates = NULL, calibration_args = list() )
models |
A named list of prediction objects, such as those returned
by
For external predictions, construct this format as shown in the example.
The separate |
metrics |
Character vector of metrics to compute (passed through to
|
t_star |
(Optional) Evaluation time for Brier score, AUC, point R2, and calibration. Defaults to the median observed follow-up. |
tau |
(Optional) The maximum prediction time horizon for R2 and C index. Default is maximum of observation time. |
event_type |
Integer code for the cause of interest (passed through). |
digits |
Number of decimal places to round the metric values (default = 2). |
censoring |
Censoring estimator for generalized R-squared metrics,
|
censor_covariates |
Optional numeric censoring-model covariates shared
across models. Rows must match each model's subject order. Passed to
|
calibration_args |
Optional named list passed to |
See pa_survival_eval_cr() for the event-observation weighting
rule, original-time support requirements and probability floor.
A data.frame with:
each row = metric
each column = model name
cell values = metric values
# Competing risks: cause 1 is the event of interest, cause 2 the competing # risk. Both cause-specific hazards are modelled. cr <- pa_simulate_fine_gray(n = 150, v = 2, beta1 = c(0.5, -0.5), beta2 = c(-0.3, 0.3), censor = 0.3, seed = 2025) cr$time <- cr$obs.times cr$status <- cr$obs.event m1 <- survival::coxph(survival::Surv(time, status == 1) ~ X1 + X2, data = cr, x = TRUE, y = TRUE) m2 <- survival::coxph(survival::Surv(time, status == 2) ~ X1 + X2, data = cr, x = TRUE, y = TRUE) p <- pa_predict_cif(model1 = m1, model2 = m2, newdata = cr, covs = c("X1", "X2"), event.type = 1, tau = max(cr$time)) pa_summarize_cr(list(csh = p), event_type = 1) # External CIF matrix: rows are grid times, columns are subjects. external_cif <- p$cif_pred[, -1, drop = FALSE] external_grid <- p$cif_pred[, 1] external <- list(cif_pred = cbind(time = external_grid, external_cif), times = p$times, status = p$status) pa_summarize_cr(list(external = external), metrics = "brier_score")# Competing risks: cause 1 is the event of interest, cause 2 the competing # risk. Both cause-specific hazards are modelled. cr <- pa_simulate_fine_gray(n = 150, v = 2, beta1 = c(0.5, -0.5), beta2 = c(-0.3, 0.3), censor = 0.3, seed = 2025) cr$time <- cr$obs.times cr$status <- cr$obs.event m1 <- survival::coxph(survival::Surv(time, status == 1) ~ X1 + X2, data = cr, x = TRUE, y = TRUE) m2 <- survival::coxph(survival::Surv(time, status == 2) ~ X1 + X2, data = cr, x = TRUE, y = TRUE) p <- pa_predict_cif(model1 = m1, model2 = m2, newdata = cr, covs = c("X1", "X2"), event.type = 1, tau = max(cr$time)) pa_summarize_cr(list(csh = p), event_type = 1) # External CIF matrix: rows are grid times, columns are subjects. external_cif <- p$cif_pred[, -1, drop = FALSE] external_grid <- p$cif_pred[, 1] external <- list(cif_pred = cbind(time = external_grid, external_cif), times = p$times, status = p$status) pa_summarize_cr(list(external = external), metrics = "brier_score")
This function calculates various predictive performance metrics for survival models based on predicted survival probabilities and observed survival data. It supports a range of evaluation measures, including explained variation, concordance indices, Brier Score, and time-dependent AUC.
pa_survival_eval( model, event_time, predicted_probability, pred_mean_survival = NULL, status, covariates, new_data = NULL, metrics = NULL, t_star = NULL, tau = NULL, calibration_args = list() )pa_survival_eval( model, event_time, predicted_probability, pred_mean_survival = NULL, status, covariates, new_data = NULL, metrics = NULL, t_star = NULL, tau = NULL, calibration_args = list() )
model |
A fitted survival model object. Optional for the metrics
computed directly from |
event_time |
A numeric vector of observed survival times. |
predicted_probability |
Numeric matrix of predicted survival probabilities, with
|
pred_mean_survival |
Optional numeric vector of predicted mean survival times, used for explained variation and concordance. Predicted probabilities are still needed for point metrics, Brier score, AUC, and calibration. |
status |
A numeric vector indicating event occurrence (1 for event, 0 for censoring). |
covariates |
A character vector specifying the names of the covariates used in the model. |
new_data |
Optional data frame used by metrics that require refitting or
prediction from |
metrics |
A character vector specifying the evaluation metrics to compute. Options include:
The default |
t_star |
(Optional) Evaluation time for point R-squared, Brier score, time-dependent AUC, and calibration. If NULL, it defaults to the median observed follow-up time. |
tau |
(Optional) A numeric value specifying the truncation time for calculating explained variation metrics. If NULL, it defaults to the maximum observed survival time. |
calibration_args |
Optional named list passed to |
Point R-squared uses the event indicator I(T <= t_star) and
the predicted event probability 1 - S at the selected prediction-grid
column. Events at the evaluation time are included.
Calibration uses 1 - S(t_star) from the last supplied grid column
at or before t_star, interpreting the grid as a right-continuous
step function (risk zero before its first time). Its outcome horizon is
exactly t_star, not the nearest future grid time. See
pa_calibration() for the estimator and boundary-probability policy.
A failed calibration fit warns and returns missing calibration metrics;
other requested metrics are retained.
A data frame with two columns:
Metric |
The name of the computed evaluation metric. |
Value |
The corresponding computed value. |
Schemper, M. and R. Henderson (2000). Predictive accuracy and explained variation in Cox regression. Biometrics 56, 249–255.
Lusa, L., R. Miceli and L. Mariani (2007). Estimation of predictive accuracy in survival analysis using R and S-PLUS. Computer Methods and Programs in Biomedicine 87, 132–137.
Potapov, S., Adler, W., Schmid, M., Bertrand, F. (2024). survAUC: Estimating Time-Dependent AUC for Censored Survival Data. R package version 1.3-0. DOI: doi:10.32614/CRAN.package.survAUC. Available at https://CRAN.R-project.org/package=survAUC.
# A Cox model evaluated on right-censored data. The simulator is seeded, # so the reported metric values are reproducible. d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # the default metric set pa_survival_eval( model = fit, event_time = pred$times, predicted_probability = pred$surv_prob, status = pred$status, covariates = c("x1", "x2"), new_data = d ) # a chosen subset; see pa_metric_names() for the accepted identifiers pa_survival_eval( model = fit, event_time = pred$times, predicted_probability = pred$surv_prob, status = pred$status, covariates = c("x1", "x2"), new_data = d, metrics = c("harrell_c", "brier_score") )# A Cox model evaluated on right-censored data. The simulator is seeded, # so the reported metric values are reproducible. d <- pa_sim_cox_weibull(n = 150, pi_c = 0.3, v = 2, beta = c(0.5, -0.5), seed = 2025) d <- d[, c("time", "status", "x1", "x2")] fit <- survival::coxph(survival::Surv(time, status) ~ x1 + x2, data = d, x = TRUE, y = TRUE) pred <- pa_predict_coxph(model = fit, covs = c("x1", "x2"), new_data = d) # the default metric set pa_survival_eval( model = fit, event_time = pred$times, predicted_probability = pred$surv_prob, status = pred$status, covariates = c("x1", "x2"), new_data = d ) # a chosen subset; see pa_metric_names() for the accepted identifiers pa_survival_eval( model = fit, event_time = pred$times, predicted_probability = pred$surv_prob, status = pred$status, covariates = c("x1", "x2"), new_data = d, metrics = c("harrell_c", "brier_score") )
This function evaluates the predictive performance of a survival model in the presence of competing risks. It calculates various performance metrics including Pseudo R-squared, R-squared, L-squared, C-index, Brier Score, and Time-Dependent AUC.
pa_survival_eval_cr( pred_cif, event_time, time.cif, status, metrics = NULL, t_star = NULL, tau = NULL, event_type = 1, censoring = c("cox", "km"), censor_covariates = NULL, calibration_args = list() )pa_survival_eval_cr( pred_cif, event_time, time.cif, status, metrics = NULL, t_star = NULL, tau = NULL, event_type = 1, censoring = c("cox", "km"), censor_covariates = NULL, calibration_args = list() )
pred_cif |
Numeric matrix of cumulative incidence functions with one
row per value of |
event_time |
A numeric vector of observed event or censoring times for subjects. |
time.cif |
A numeric vector of time points corresponding to the predicted CIFs. |
status |
A numeric vector indicating the event type for each subject: by default
|
metrics |
A character vector specifying the performance metrics to calculate.
The default
|
t_star |
Optional numeric value specifying a reference time for evaluating point-version Pseudo_R_square,
the Brier score, time-dependent AUC, and calibration. Defaults to the median
of |
tau |
Optional numeric value specifying a restricted time horizon for model evaluation. Defaults to the maximum observed follow-up time. |
event_type |
An integer indicating the primary event type of interest (default = |
censoring |
Censoring estimator for the generalized R-squared metrics:
|
censor_covariates |
Optional numeric vector, matrix or data frame of
censoring-model covariates, in the same subject order as |
calibration_args |
Optional named list passed to |
The generalized R-squared metrics use event-observation IPCW. Censoring is
estimated on the original follow-up times, treating any observed event type
as an observed event. The unnormalized weight is
; weights are normalized to sum
to one. A censored subject has zero weight even if follow-up exceeds tau.
Restricting the working outcome therefore does not limit censoring support
requirements to tau: support at contributing original event times is needed.
Cox weights require conditional independent censoring and a valid censoring
model; marginal weights require the corresponding marginal independence.
Censoring survival probabilities are floored at .
The interval outcome is the target-event time capped at tau, with competing
events assigned to tau. The point outcome is I(T <= t_star, D = event_type),
including target events at the evaluation time, as in the CIF. Normalized ratio
estimates are not asserted to be exactly unbiased in finite samples.
Calibration uses the target-event CIF at or before t_star on the
supplied grid and evaluates the outcome at t_star. It uses the local
Aalen–Johansen estimator described in pa_calibration_cr(); censoring
and censor_covariates do not control this estimator. A failed curve
fit warns and returns missing calibration metrics without discarding others.
A data frame with two columns:
"Metric": Name of the evaluated metric.
"Value": Computed value of the corresponding metric.
# Competing risks: cause 1 is the event of interest, cause 2 the competing # risk. Both cause-specific hazards are modelled. cr <- pa_simulate_fine_gray(n = 150, v = 2, beta1 = c(0.5, -0.5), beta2 = c(-0.3, 0.3), censor = 0.3, seed = 2025) cr$time <- cr$obs.times cr$status <- cr$obs.event m1 <- survival::coxph(survival::Surv(time, status == 1) ~ X1 + X2, data = cr, x = TRUE, y = TRUE) m2 <- survival::coxph(survival::Surv(time, status == 2) ~ X1 + X2, data = cr, x = TRUE, y = TRUE) p <- pa_predict_cif(model1 = m1, model2 = m2, newdata = cr, covs = c("X1", "X2"), event.type = 1, tau = max(cr$time)) pa_survival_eval_cr( pred_cif = p$cif_pred[, -1], event_time = p$times, time.cif = p$cif_pred[, 1], status = p$status, event_type = 1 )# Competing risks: cause 1 is the event of interest, cause 2 the competing # risk. Both cause-specific hazards are modelled. cr <- pa_simulate_fine_gray(n = 150, v = 2, beta1 = c(0.5, -0.5), beta2 = c(-0.3, 0.3), censor = 0.3, seed = 2025) cr$time <- cr$obs.times cr$status <- cr$obs.event m1 <- survival::coxph(survival::Surv(time, status == 1) ~ X1 + X2, data = cr, x = TRUE, y = TRUE) m2 <- survival::coxph(survival::Surv(time, status == 2) ~ X1 + X2, data = cr, x = TRUE, y = TRUE) p <- pa_predict_cif(model1 = m1, model2 = m2, newdata = cr, covs = c("X1", "X2"), event.type = 1, tau = max(cr$time)) pa_survival_eval_cr( pred_cif = p$cif_pred[, -1], event_time = p$times, time.cif = p$cif_pred[, 1], status = p$status, event_type = 1 )