At first, I would like to plot the survival curves. After that , the main
use will be to calculate conditional probabilities - given that an
individual already survived x days, what will be the chance it survives
till day ,e.g., 100.
I 've already found a solution, though.
Of the - non-subscribed - survfit-object, the following functions select
the indices belonging to the individual.
Usage is e.g.:
stratas= extract_strata(coxph.3, data)
coxph.3.surfvit = survfit(coxph.3, newdata=data)
strata_subscripts = extract_strata_subscripts(coxph.3.surfvit)
plot(0,0, ylim=c(0,1), xlim=c(-10,360))
for(i in 1:100) lines(coxph.3.surfvit$time[strata_subscripts[,stratas[i]]
],
coxph.3.surfvit$surv[strata_subscripts[,stratas[i]]
,i],
col = rainbow(100)[i])
Of course, this is not very effective, but it works. Optimisations could
be to work with basehaz and predict, instead of fitting a survfit-object.
But, here are the functions I use (still partly in german, nor well
documented, and no nice coding)
All the best
Julian
extract_strata <- function( object, data){
## returns the correct stratum for every element of data
## Gibt für jedes Element von Data das zugehörige Stratum gemäß object
zurück
if( inherits(object,"coxph")){
terms_form = terms(object$formula, specials="strata")
} else if (inherits(object,"formula")){
terms_form = terms(object, specials="strata")
} else stop("object muss von Klasse coxph oder formula sein")
if(length(attr(terms_form, which="specials")$strata)==0) {
warning("Keine strata gefunden")
return (NULL)
} else if(length(attr(terms_form, which="specials")$strata)>1) {
warning("Mehr als ein Aufruf von strata gefunden, return NULL")
return (NULL)
} else {
# strata-call parsen, im envrionment data ausführen, zurückgeben)
return(eval(parse(text=rownames(
attr(terms_form,"factors"))[attr(terms_form, which="specials")$strata]),
envir=data) )
}
}
extract_strata_subscripts <- function(survfit_object){
if( !inherits(survfit_object,"survfit.cox")) stop("survfit_object must
be of class \"survfit.cox\"")
require(survival)
if( is.null(survfit_object$strata)){
warning("No Strata found")
return(TRUE)
}
stratanames= names(survfit_object$strata)
nstrata = length(stratanames)
ntimes = length(survfit_object$time)
strataborders = matrix(ncol=3, nrow=nstrata,
dimnames=list(
strata=stratanames,borders=c("min","max", "length")))
strataborders[1,]=c(1,nrow(survfit_object[1]$surv),nrow(survfit_object[1]$
surv))
for(x in 2:nstrata){
strataborders[x,1]=strataborders[x-1,2]+1
strataborders[x,2]=strataborders[x-1,2]+ nrow(survfit_object[x]$surv)
strataborders[x,3] = nrow(survfit_object[x]$surv)
}
ret_matrix=matrix(data=F,ncol=nstrata, nrow=ntimes)
colnames(ret_matrix)=stratanames
attr(ret_matrix, which="strataborders")= strataborders
for(x in stratanames)
ret_matrix[( (1:ntimes)>= strataborders[x,1])&( (1:ntimes)<=
strataborders[x,2]),x] =T
return(ret_matrix)
}
-----Ursprüngliche Nachricht-----
Von: Terry Therneau [mailto:thern...@mayo.edu]
Gesendet: Freitag, 26. Juli 2013 16:12
An: r-help@r-project.org; julian.bo...@elitepartner.de
Betreff: Re: [R] prediction survival curves for coxph-models; how to
extract the right strata per individual
It would help me to give advice if I knew what you wanted to do with the
new curves.
Plot, print, extract?
A more direct solution to your question will appear in the next release of
the code, btw.
Terry T.
On 07/25/2013 05:00 AM, r-help-requ...@r-project.org wrote:
> My problem is:
>
>
>
> I have a coxph.model with several strata and other covariables.
>
> Now I want to fit the estimated survival-curves for new data, using
> survfit.coxph.
>
> But this returns a prediction for each stratum per individual. So if I
> have 15 new individuals and 10 strata, I have 150 fitted surivival
> curves (or if I don't use the subscripts I have 15 predictions with
> the curves for all strata pasted together)
>
>
>
> Is there any possibility to get only the survival curves for the
> stratum the new individual belongs to?
>
>
______________________________________________
R-help@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.
