[R] Obtaining value of median survival for survfit function to use in calculation

[Frank E Harrell Jr]

Or try the following:

require(rms)
data(lung)
survfit(Surv (time, status) ~ ph.ecog, data=lung)
lung <- transform(lung, ph.ecog=factor(ph.ecog))
f <- cph(Surv(time, status) ~ strat(ph.ecog), data=lung, surv=TRUE)
## doesn't assume PH
quant <- Quantile(f)   # print quant to see how all quantiles are computed
quant(stratum=1)
meds <- sapply(1:4, function(x) quant(stratum=x))
meds

0.5 0.5 0.5 0.5
394 306 199 118

## Much more precise, if PH holds:
dd <- datadist(lung); options(datadist='dd')
f <- cph(Surv(time, status) ~ ph.ecog, data=lung, surv=TRUE)
quant <- Quantile(f)
Predict(f, ph.ecog=., fun=function(lp) quant(lp=lp))
## Note that confidence intervals are approximate because they don't
## take into account uncertainty in the underlying survival curve

   ph.ecog yhat lower upper
1       0  429   583   320
2       1  329   329   329
3       2  210   288   176
4       3   88   353    11

## Test for linearity in ph.ecog
f <- cph(Surv(time, status) ~ scored(ph.ecog), data=lung)
anova(f)  # chi-sq=1.05 with 2 d.f. P=.59
## Note this type of testing distorts the final d.f.

                 Wald Statistics          Response: Surv(time, status)

  Factor     Chi-Square d.f. P
  ph.ecog    20.27      3    0.0001
   Nonlinear  1.05      2    0.5925
  TOTAL      20.27      3    0.0001

## More precise, as saves two d.f. by assuming linearity, PH
lung <- transform(lung, ph.ecog.num=as.numeric(as.character(ph.ecog)))
dd <- datadist(lung)
f <- cph(Surv(time, status) ~ ph.ecog.num, data=lung, surv=TRUE)
quant <- Quantile(f)
Predict(f, ph.ecog.num=., fun=function(lp) quant(lp=lp))
## Same caution about confidence limits

   ph.ecog.num yhat lower upper
1           0  450   550   364
2           1  310   310   310
3           2  212   269   183
4           3  166   210   122

quant(lp=predict(f, data.frame(ph.ecog.num=0:3)))

   1   2   3   4
450 310 212 166

Frank


David Winsemius wrote:
> The survfit.object help page says:
> 
> "The print.survfit method does more computation than is typical for a 
> print method and is documented on a separate page."
> 
> It takes a bit of digging, but after first trying:
> 
> getAnywhere(print.survfit)   # and then following code and trying
> getAnywhere(survmean)  # survmean is the function which does the work 
> using pfun
> 
> ... I think I finally understand exactly what the print.survfit help 
> page means when it refers to a "side-effect".
> 
> So there is no durable element in the survfit object that is the median. 
> The automation you request would involve duplicating the code of 
> survmean but with assignment of the "out" matrix to something that does 
> not get discarded.
> 


-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University