riskRegression (version 2025.09.17)

predict.CauseSpecificCox: Predicting Absolute Risk from Cause-Specific Cox Models

Description

Apply formula to combine two or more Cox models into absolute risk (cumulative incidence function).

Usage

# S3 method for CauseSpecificCox
predict(
  object,
  newdata,
  times,
  cause,
  type = "absRisk",
  landmark = NA,
  keep.times = 1L,
  keep.newdata = 1L,
  keep.strata = 1L,
  se = FALSE,
  band = FALSE,
  iid = FALSE,
  confint = (se + band) > 0,
  average.iid = FALSE,
  product.limit = TRUE,
  store = NULL,
  diag = FALSE,
  max.time = NULL,
  ...
)

Arguments

object

The fitted cause specific Cox model

newdata

[data.frame or data.table] Contain the values of the predictor variables defining subject specific predictions relative to each cause. Should have the same structure as the data set used to fit the object.

times

[numeric vector] Time points at which to return the estimated absolute risk.

cause

[integer/character] Identifies the cause of interest among the competing events.

type

[character] Can be changed to "survival" if the event free survival should be output instead of the absolute risk.

landmark

[integer] The starting time for the computation of the cumulative risk.

keep.times

[logical] If TRUE add the evaluation times to the output.

keep.newdata

[logical] If TRUE add the value of the covariates used to make the prediction in the output list.

keep.strata

[logical] If TRUE add the value of the strata used to make the prediction in the output list.

se

[logical] If TRUE compute and add the standard errors to the output.

band

[logical] If TRUE compute and add the quantiles for the confidence bands to the output.

iid

[logical] If TRUE compute and add the influence function to the output.

confint

[logical] If TRUE compute and add the confidence intervals/bands to the output. They are computed applying the confint function to the output.

average.iid

[logical]. If TRUE add the average of the influence function over newdata to the output.

product.limit

[logical]. If TRUE the survival is computed using the product limit estimator. Otherwise the exponential approximation is used (i.e. exp(-cumulative hazard)).

store

[vector of length 2] Whether prediction should only be computed for unique covariate sets and mapped back to the original dataset (data="minimal") and whether the influence function should be stored in a memory efficient way (iid="minimal"). Otherwise use data="full" or iid="full".

diag

[logical] when FALSE the absolute risk/survival for all observations at all times is computed, otherwise it is only computed for the i-th observation at the i-th time.

max.time

[numeric] maximum time of the response of the fitted data. Only relevant if model response element has been removed

...

not used.

Author

Brice Ozenne [e-mail], Thomas A. Gerds [e-mail]

Details

This function computes the absolute risk as given by formula 2 of (Ozenne et al., 2017). Confidence intervals and confidence bands can be computed using a first order von Mises expansion. See the section "Construction of the confidence intervals" in (Ozenne et al., 2017).

A detailed explanation about the meaning of the argument store = c(iid="full") can be found in (Ozenne et al., 2017) Appendix B "Saving the influence functions".

The function is not compatible with time varying predictor variables nor frailty.

The iid decomposition is output using an array containing the value of the influence of each subject used to fit the object (dim 1), for each subject in newdata (dim 3), and each time (dim 2).

References

Brice Ozenne, Anne Lyngholm Sorensen, Thomas Scheike, Christian Torp-Pedersen and Thomas Alexander Gerds. riskRegression: Predicting the Risk of an Event using Cox Regression Models. The R Journal (2017) 9:2, pages 440-460.

See Also

confint.predictCSC to compute confidence intervals/bands. autoplot.predictCSC to display the predictions.

Examples

library(survival)
library(prodlim)
library(data.table)

#######################
#### generate data ####
#######################

set.seed(5)
d <- sampleData(80,outcome="comp") ## training dataset
nd <- sampleData(4,outcome="comp") ## validation dataset
d$time.round <- round(d$time,1) ## create tied events
ttt <- sort(sample(x = unique(d$time), size = 10))

########################
#### Aalen-Johansen ####
########################

#### no ties
fit.AE <- CSC(Hist(time, event) ~ 1, data = d)
GS.AE <- prodlim(Hist(time, event) ~ 1, data = d)

## product limit
ePL.AE <- predict(fit.AE, newdata = d[1], times = GS.AE$time,
                  product.limit = TRUE, cause = 1)
range(ePL.AE$absRisk - GS.AE$cuminc[[1]]) ## same

## exponential approximation
predict(fit.AE, newdata = d[1], times = GS.AE$time,
        product.limit = FALSE, cause = 1)

#### with ties (use Breslow estimator instead of Efron to match prodlim)
fit.AE.ties <- CSC(Hist(time.round, event) ~ 1, ties = "breslow",
                   data = as.data.frame(d)[c("time.round","event")])
GS.AE.ties <- prodlim(Hist(time.round, event) ~ 1,
                   data = as.data.frame(d)[c("time.round","event")])

## product limit
ePL.AE.ties <- predict(fit.AE.ties, newdata = d[1], times = GS.AE.ties$time,
                       product.limit = TRUE, cause = 1)
range(ePL.AE.ties$absRisk - GS.AE.ties$cuminc[[1]]) ## same

###################################
#### Stratified Aalen-Johansen ####
###################################

#### no ties
fit.SAE <- CSC(Hist(time, event) ~ strata(X1), data = d)
GS.SAE <- prodlim(Hist(time, event) ~ X1, data = d)

## product limit
index.strata1 <- GS.SAE$first.strata[1]:(GS.SAE$first.strata[2]-1)
ePL.SAE <- predict(fit.SAE, newdata = d[1],
                   times = GS.SAE$time[index.strata1],
                   product.limit = TRUE, cause = 1)
range(ePL.SAE$absRisk - GS.SAE$cuminc[[1]][index.strata1]) ## same

## exponential approximation
predict(fit.SAE, newdata = d[1:10], times = 1:10, product.limit = FALSE)

############################
#### Cause-specific Cox ####
############################

## estimate a CSC model based on the coxph function
CSC.fit <- CSC(Hist(time,event)~ X3+X8, data=d, method = "breslow")

## compute the absolute risk of cause 1, in the validation dataset
## at time 1:10
CSC.risk <-  predict(CSC.fit, newdata=nd, times=1:10, cause=1)
CSC.risk

## compute absolute risks with CI for cause 2
## (without displaying the value of the covariates)
predict(CSC.fit,newdata=nd,times=1:10,cause=2,se=TRUE,
        keep.newdata = FALSE)

## other example
library(survival)
CSC.fit.s <- CSC(list(Hist(time,event)~ strata(X1)+X2+X9,
 Hist(time,event)~ X2+strata(X4)+X8+X7),data=d, method = "breslow")
predict(CSC.fit.s,cause=1,times=ttt,se=1L) ## note: absRisk>1 due to small number of observations

## using the cph function instead of coxph
CSC.cph <- CSC(Hist(time,event)~ X1+X2,data=d, method = "breslow", fitter = "cph")#' 
predict(CSC.cph, newdata = d, cause = 2, times = ttt)

## landmark analysis
T0 <- 1
predCSC.afterT0 <- predict(CSC.fit, newdata = d, cause = 2, times = ttt[ttt>T0], landmark = T0)
predCSC.afterT0