[R] prediction survival curves for coxph-models; how to extract the right strata per individual

[julian.bothe at elitepartner.de]

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:therneau at mayo.edu]
Gesendet: Freitag, 26. Juli 2013 16:12
An: r-help at r-project.org; julian.bothe at 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-request at 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?
>
>