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
______________________________________________
[email protected] 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.
