[R] Hazard rate estimation by smoothing baseline cumulative hazard from Cox model - Was: RE: Competing risks Kalbfleisch & Prentice method
Ravi Varadhan
RVaradhan at jhmi.edu
Thu Mar 26 21:26:20 CET 2009
Hi Eleni,
I will take a look at this. I have some preliminary comments.
You estimate the hazard function from the Cox model baseline cumulative
hazard by differencing successive jumps. It seems that a better approach
might be to estimate this using kernel smoothing, i.e. as the derivative of
kernel-smoothed cumulative hazard function. This method is available in the
"muhaz" package. However, the muhaz() function does not work with the Cox
model baseline cumulative hazard. It requires you to input the original
data on times and censoring indicators. It would be nice if this were
possible. I am wondering why Terry Therneau's "survival" package doesn't
have this option.
Best,
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of Eleni Rapsomaniki
Sent: Thursday, March 26, 2009 1:17 PM
To: Ravi Varadhan; Heinz Tuechler
Cc: r-help at r-project.org
Subject: Re: [R] Competing risks Kalbfleisch & Prentice method
Ravi,
I agree with you, that the Fine & Gray method does much more than
calculating the cumulative incidence. The Kalbfleisch & Prentice formula
relies on the strong assumption that the coefficients for the cause specific
hazard models are valid in the presence of competing risks. You will find my
code below, but it may be wrong (I'm not really a statistician!). And no, it
doesn't even calculate confidence intervals...
Thank you for pointing out these papers.
Please let me know (nicely) if you find bugs!!!
#####
#Method to calculate Cumulative Incidence adjusting for competing risks,
based on the Kalbfleisch & Prentice formula, 1980, p.169 #create data ftime
<- rexp(200, 0.3) fstatus <- sample(0:2,200,replace=TRUE) cov <-
matrix(runif(600),nrow=200) dimnames(cov)[[2]] <- c('x1','x2','x3')
dat=data.frame(ftime, fstatus, cov)
diffrHaz=function(x){
hazard=x$hazard
time=x$time
hz.v=NULL
for(e in 1:(length(hazard)-1)){
hz=hazard[e+1]-hazard[e]
hz.v=c(hz.v,hz)
}
dhz=data.frame(hazard=hz.v, time=x$time[-length(x$time)])
return(dhz)
}
#build cause specific cox models
library(Design)
a.f=cph(Surv(ftime, fstatus==1) ~ x1+x2+x3, data=dat, surv=T, x=T, y=T)
b.f=cph(Surv(ftime, fstatus==2) ~ x1+x2+x3, data=dat, surv=T, x=T, y=T)
#get unique event times - up till the time of interest, eg. time=10
uts=unique(dat$ftime[dat$fstatus>0] ) uts=uts[uts<=10] # the times up till
10 years uts=uts[order(uts)] #order them
#get baseline hazard rate (cumulative)
a.bz=basehaz(a.f)
b.bz=basehaz(b.f)
#get instanteneous baseline hazard rate, apply above function
a.dhz1=diffrHaz(a.bz)
b.dhz1=diffrHaz(b.bz)
P.t=NULL #intitialize where to store all results for(i in 1:nrow(dat)){ # i
is the individual, iterate through the number of rows in dat
P.i.t=NULL
all.lambda.i.t=NULL
for(u.i in 1:length(uts)){
u=uts[u.i]
#get the instanteneous baseline hazard from each model corresponding
to this timepoint (if there is no event at that time return 0)
a.dhz.u=ifelse((u %in%
a.dhz1$time)==F,0,a.dhz1$hazard[a.dhz1$time==u])
b.dhz.u=ifelse((u %in%
b.dhz1$time)==F,0,b.dhz1$hazard[b.dhz1$time==u])
#multiply by the linear predictors to get the actual hazard at that
point
a.lambda.i.u=a.dhz.u*exp(a.f$linear.predictors[i])
b.lambda.i.u=b.dhz.u*exp(b.f$linear.predictors[i])
all.lambda.i.u=a.lambda.i.u+b.lambda.i.u
#store in a vector
all.lambda.i.t=c(all.lambda.i.t, all.lambda.i.u) #build a vector
with all previous all.lambda.i.u's
S.i.u=exp(-sum(all.lambda.i.t)) #this is the probability of
surviving till time u
p.i.u=S.i.u*a.lambda.i.u
P.i.t=c(P.i.t, p.i.u)
}
P.t[i]=sum(P.i.t)
}
# to compare with unadjusted:
a.risk = 1-survest(a.f, linear.predictors=a.f$linear.predictors,
times=10)$surv
plot(a.risk ~ P.t) #P.t is the CR adjusted, should always be lower
mean(a.risk)
mean(P.t)
Eleni Rapsomaniki
Research Associate
Tel: +44 (0) 1223 740273
Strangeways Research Laboratory
Department of Public Health and Primary Care University of Cambridge
-----Original Message-----
From: Ravi Varadhan [mailto:RVaradhan at jhmi.edu]
Sent: 26 March 2009 14:36
To: Eleni Rapsomaniki; 'Arthur Allignol'
Cc: r-help at r-project.org
Subject: RE: [R] Competing risks Kalbfleisch & Prentice method
Hi Eleni,
I would like to take a look at your R function for obtaining the cumulative
incidence function (CIF) from individual Cox models for cause-specific
hazards (CSH). Does your code predict the CIF (with pointwise confidence
intervals and global confidence bands) for different sub-groups? Have you
seen the paper by Cheng, Fine, and Wei (Biometrics 1998) that does this?
A major advantage of the F&G model is that you can get a direct, numerical
measure of the effect of a covariate on the CIF. This cannot be obtained by
modeling all the CSHs and then combining them. The idiosyncratic assumption
concerning risk set in F&G model is made mainly for mathematical purposes so
that a proportional hazards form may be obtained for the CIF. You can test
this assumption by plotting schonefeld-type residuals (this is available in
cmprsk). Fine (Biostatistics 2006) provides a different approach that
relaxes this assumption (it also uses a different estimation approach), but
I don't know if there is an R implementation for that.
Thanks,
Ravi.
------------------------------------------------------------------------
----
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
------------------------------------------------------------------------
----
--------
-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
Behalf Of Eleni Rapsomaniki
Sent: Thursday, March 26, 2009 10:18 AM
To: Arthur Allignol
Cc: r-help at r-project.org
Subject: Re: [R] Competing risks Kalbfleisch & Prentice method
Thank you for your reply.
It wasn't too hard to code actually, which is probably why it doesn't have a
special package dedicated to it. The results are almost identical to Fine &
Gray regression model. The problem with the latter is that my colleagues are
not convinced that the model assumptions (people who die from competing
causes remaining in the risk set) are theoretically sound.
If anybody is interested in the Kalbfleisch & Prentice based cumulative
incidence adjusting for competing risks with covariates, I'm happy to supply
the code.
Eleni Rapsomaniki
Research Associate
Tel: +44 (0) 1223 740273
Strangeways Research Laboratory
Department of Public Health and Primary Care University of Cambridge
-----Original Message-----
From: Arthur Allignol [mailto:arthur.allignol at fdm.uni-freiburg.de]
Sent: 26 March 2009 10:36
To: Eleni Rapsomaniki
Cc: r-help at r-project.org
Subject: Re: [R] Competing risks Kalbfleisch & Prentice method
I don't think there is a package to do that.
But you could have a look at ?predict.crr.
Best regards,
Arthur Allignol
Eleni Rapsomaniki wrote:
>
>
> Dear R users
>
>
>
> I would like to calculate the Cumulative incidence for an event
> adjusting for competing risks and adjusting for covariates. One way to
> do this in R is to use the cmprsk package, function crr. This uses the
> Fine & Gray regression model. However, a simpler and more classical
> approach would be to implement the Kalbfleisch & Prentice method
(1980,
> p 169), where one fits cause specific cox models for the event of
> interest and each type of competing risk, and then calculates
incidence
> based on the overall survival. I believe that this is what the cuminc
> function in the aforementioned package does, but it does not allow to
> adjust for a vector of covariates.
>
>
>
> My question is, is there an R package that implements the Kalbfleisch
&
> Prentice method for competing risks with covariates?
>
>
>
> for example, if k1 is the cause of interest among k competing causes:
>
> P_k1(t; x)=P(T<=t, cause=k1|x)=Sum(u=0, ..., u=t)
{hazard_k(u;x)*S(u;x)}
>
> where S(u;x) = exp{-sum_of_k(sum(hazard_k(u))}
>
>
>
> I have searched extensively for an implementation of this in many
> packages, but it appears that more complex approaches are more
commonly
> implemented, such as timereg package.
>
>
>
> Eleni Rapsomaniki
>
>
>
> Research Associate
>
> Strangeways Research Laboratory
>
> Department of Public Health and Primary Care
>
>
>
> University of Cambridge
>
>
>
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
More information about the R-help
mailing list