I'd like to calculate the expected risk (cumulative incidences), which are derived from a 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 the expected risk (cumulative incidence of events) for each case according to X1 = 0 or X1 = 1.
In other words, if the X1 of all 10,000 cases were 0, how could I calculate the expected risk for each case? At the same time, if the X1 of all 10,000 cases were 1, how could I calculate the expected risk for each case? (using R)
After calculating the expected risk for each patient, then I'd like to calculate the risk-benefit ratio for each patient according to variable X1 [by (Expected risk, when X1 = 0) - (Expected risk, when X1 = 1)], does 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.