[R] Calculating the probability of an event at time "t" from a Cox model fit
David Winsemius
dwinsemius at comcast.net
Mon Dec 19 16:11:08 CET 2011
On Dec 19, 2011, at 3:34 AM, aajit75 wrote:
> Dear R-users,
>
> I would like to determine the probability of event at specific time
> using
> cox model fit. On the development sample data I am able to get the
> probability of a event at time point(t).
> I need probability score of a event at specific time, using scoring
> scoring
> dataset which will have only covariates and not the response
> variables.
>
> Here is the sample code:
>
> n = 1000
> beta1 = 2; beta2 = -1;
> lambdaT = .02 # baseline hazard
> lambdaC = .4 # hazard of censoring
> x1 = rnorm(n,0)
> x2 = rnorm(n,0)
>
> # true event time
> T = rweibull(n, shape=1, scale=lambdaT*exp(-beta1*x1-beta2*x2))
> C = rweibull(n, shape=1, scale=lambdaC) #censoring time
> time = pmin(T,C) #observed time is min of censored and true
> event = time==T # set to 1 if event is observed
> dataphr=data.frame(time,event,x1,x2)
>
> library(survival)
> fit_coxph <- coxph(Surv(time, event)~ x1 + x2 , method="breslow")
>
> library(peperr)
> predictProb.coxph(fit_coxph, Surv(dataphr$time, dataphr$event),
> dataphr,
> 0.003)
>
> # Using predictProb.coxph function, probability of event at time (t)
> is
> estimated for cox fit models,
No , it should not be the probability at time T (which would be
infinitesimally small at any time point), but rather the probability
prior to time T. I am puzzled why one would not use predict.coxph with
type="expected".
?predict.coxph
> I want to estimate this probability on scoring
> dataset score_data as below with covariate x1 and x2.
>
> Is it possible/ is there any way to get these probabilities? since in
> predictProb.coxph function it requires response, which is not
> preseent on
> scoring sample.
Raising what would seem to be the obvious question: Why didn't you
make one? Or use the derivation Surv object?
>
> n = 1000
> set.seed(1)
> x1 = rnorm(n,0)
> x2 = rnorm(n,0)
> score_data <- data.frame(x1,x2)
> n = 100
> set.seed(1)
> x1 = rnorm(n,0)
> x2 = rnorm(n,0)
> score_data <- data.frame(x1,x2)
> s.scr <- survival::Surv(rep(2, n),rep(1, n) )
> predictProb.coxph(fit_coxph, s.scr, score_data, 0.003)
Error in model.frame.default(formula = Surv(time, event) ~ x1 + x2) :
variable lengths differ (found for 'x1')
After getting that model.frame error with the 10000 element test set:
"variable lengths differ (found for 'x1')", I changed it to the same
length as in the derivation set and ran:
predictProb.coxph(fit_coxph, Surv(time, event), score_data, 0.003)
[,1]
[1,] 9.977413e-01
[2,] 9.879081e-01
[3,] 9.893564e-01
[4,] 5.857123e-01
[5,] 9.550606e-01
[6,] 9.761398e-01
[7,] 9.699297e-01
snipped
So apparently that function will accept a dataframe even though its
help page says it should be a matrix and it will run if the same
number of records aare present as in the derivation set. This next
version also worked and gave identical results, so it appears the Surv
object is really just a place holder.
s.scr <- survival::Surv(rep(2, n),rep(0, n) )
> predictProb.coxph(fit_coxph, s.scr, score_data, 0.003)
[,1]
[1,] 9.977413e-01
[2,] 9.879081e-01
[3,] 9.893564e-01
[4,] 5.857123e-01
[5,] 9.550606e-01
snipped
Since I am not a regular user I cannot really answer the question
about whether a Surv() object and a dataframe of different dimensions
than the fit object should be valid input to that function.
(Copied the peperr author.)
>
David Winsemius, MD
West Hartford, CT
More information about the R-help
mailing list