How to calculate expected risk from fitted Cox PH model in R?

Question

I'd like to calculate expected risk (cumulative incidences), which are derived from fitted Cox PH model using R packages.

I have the fitted Cox PH model like as follows:

[Variables] Dataset: 10,000 cases(patients) with 6 variabels

t: time to event, s: event (coded as 0 or 1)

covariates: X1, X2, X3, X4 (coded as 0 or 1)

fit <- coxph (Surv (t,s) ~ X1+X2+X3+X4, data=data)

summary (fit)

       coef   exp(coef)  se (coef)   z     P         
X1   -0.3777   0.6855   0.1120   -3.37   0.00075
X2    0.4014   1.4938   0.0518    7.74   <0.0001
X3    0.7417   2.0995   0.0893    8.31   <0.0001
X4    0.4330   1.5419   0.1268    3.42    <0.001 

From this model, I'd like to calculate expected risk (cumulative incidence of events) for each cases according to X1 = 0 or X1 = 1.

In other words, if the X1 of all 10,000 cases were 0, how could I calculate expected risk for each cases? At the same time, if the X1 of all 10,000 cases were 1, how could I calculate expected risk for each cases? (using R)

After calculating the expected risk for each patients, then I'd like to calculate the risk-benefit ratio for each patients according to variable X1 [by (Expected risk, when X1 = 0) - (Expected risk, when X1 = 1)], dose it right?

Meanwhile, I have tried to plot the expected cumulative risk curve according to X1 like as follows; was it appropriate ?

baseha = basehaz (fit, centered=FALSE)

X1=0: exp (-0.38*0) = 1, X1=1: exp (-0.38*1) = 0.68

plot (baseha$$time[,2], baseha$hazard*(1), type="s", lty=1)

lines (baseha$$time[,2], baseha$hazard*(0.68), type="s", lty=2)

Thanks for your kind answers.

Answer 1

To derive risk estimates from a Cox model you typically use the Breslow estimator for the cumulative baseline hazard that is coupled with the hazard ratios from the Cox model. This is implemented in the survfit() function of the survival package. In particular, if you give as a first argument your fitted Cox model, and in the newdata arguments you specify a data frame with the specific values for the covariates X1, $\dots$, X4, you get survival probabilities, which you can transform to risks by calculating 1 - Survival.