[R] Nomogram with stratified cph in Design package
Frank E Harrell Jr
f.harrell at vanderbilt.edu
Sat May 2 20:33:44 CEST 2009
Renee,
I finally found the problem. You misspecified fun.at as a list with two
elements. It needed 6 elements to correspond with fun:
fun.at=list(at.surv, at.surv, at.surv,
at.med, at.med, at.med)
You have subjects with a very good prognosis and you are asking for
median survival time. The estimate of the median from the Cox model is
NA unless the linear predictor is very high. I've added some code to
understand this below and made some minor code improvements. -Frank
h <- harrell # and avoid attach
dd<- datadist(h); options(datadist="dd")
units(h$time)<-"Month"
S <- with(h, Surv(time,fail))
f<-cph(S~age+gender+size+met+node+ece+hist, x=TRUE, y=TRUE,
surv=T,time.inc=10*12,
method="breslow", data=h)
mins <- min(f$surv) # .984
sum(mins ^ exp(predict(f)) < .5) # only 241 subjects have predicted
prob < .5
srv=Survival(f)
srv120=function(lp) srv(10*12,lp)
quant=Quantile(f)
med=function(lp) quant(.5,lp)
r <- function(lp, med, srv) {
m <- med(lp)
cat(' Non-NA median estimates: ', sum(!is.na(m)),'\n',
'lp range for non-NA median: ',
range(lp[!is.na(m)]),'\n',
'range of survival prob with non-NA
median:',range(srv(lp)[!is.na(m)]),
'\n',
'median range: ', range(m,
na.rm=TRUE),'\n')
}
r(predict(f), med, srv120)
at.surv=c(0.001,0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.97,0.98,0.99,.999)
at.med=c(120,80,60,40,30,20,15,10,8,6,4,2,0)
nomogram(f, cex.var=0.75, cex.axis=0.6, ia.space=1, lp=FALSE,
fun=list(srv120, med),
funlabel=c("120-mo Survival Prob","Median Survival (months)"),
fun.at=list(at.surv, at.med),varname.label=F, maxscale=100)
title("unstratified nomogram")
f<-cph(S~age+gender+size+met+node+ece+strat(hist), x=TRUE, y=TRUE,
surv=TRUE, time.inc=10*12, method="breslow", data=h)
for(i in 1:3)
{
hi <- levels(h$hist)[i]
j <- h$hist == hi
lp <- predict(f)[j]
if(hi=='follicular') lp.f <- lp
if(hi=='other') lp.o <- lp
if(hi=='papillary') lp.p <- lp
mins <- min(f$surv[[i]])
count <- sum(mins ^ exp(lp) < .5)
prn(c(mins,count))
}
srv=Survival(f)
surv.p <- function(lp) srv(10*12, lp, stratum="hist=papillary")
surv.f <- function(lp) srv(10*12, lp, stratum="hist=follicular")
surv.o <- function(lp) srv(10*12, lp, stratum="hist=other")
quant=Quantile(f)
med.p <- function(lp) quant(.5, lp, stratum="hist=papillary")
med.f <- function(lp) quant(.5, lp, stratum="hist=follicular")
med.o <- function(lp) quant(.5, lp, stratum="hist=other")
r(lp.p, med.p, surv.p)
r(lp.f, med.f, surv.f)
r(lp.o, med.o, surv.o)
at.surv=c(0.001,0.01,0.1,0.2,0.3,0.4,0.5,06,0.7,0.8,0.9,0.95,0.97,0.98,0.99,.999)
at.med=c(120,80,60,40,30,20,15,10,8,6,4,2,0)
nomogram(f, cex.var=0.75, cex.axis=0.6, ia.space=1, lp=FALSE,
fun=list(surv.p, surv.f, surv.o, med.p, med.f, med.o),
funlabel=c("Surv Pap","Surv Fol","Surv Oth","Med surv pap",
"Med surv Fol","Med Surv Oth"),
fun.at=list(at.surv, at.surv, at.surv,
at.med, at.med, at.med), varname.label=FALSE, maxscale=100)
title("stratified nomogram")
> -----Original Message-----
> From: Frank E Harrell Jr [mailto:f.harrell at vanderbilt.edu]
> Sent: Sunday, April 26, 2009 6:13 AM
> To: David Winsemius
> Cc: reneepark; r-help at r-project.org
> Subject: Re: [R] Nomogram with stratified cph in Design package
>
> David Winsemius wrote:
> >
> > On Apr 25, 2009, at 6:57 PM, reneepark wrote:
> >
> >>
> >> Hello,
> >> I am using Dr. Harrell's design package to make a nomogram. I was
> >> able to make a beautiful one without stratifying, however, I will
> >> need to stratify to meet PH assumptions. This is where I go wrong,
> >> but I'm not sure where.
> >>
> >>
> >>
> >> Non-Stratified Nomogram:
> >>
> >> f<-cph(S~A+B+C+D+E+F+H,x=T,y=T,surv=T,time.inc=10*12,method="breslow"
> >> )
> >> srv=Survival(f)
> >> srv120=function(lp) srv(10*12,lp)
> >> quant=Quantile(f)
> >> med=function(lp) quant(.5,lp)
> >> at.surv=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
> >> at.med=c(120,80,60,40,30,20,15,10,8,6,4,2,0)
> >> nomogram(f,lp=F, fun=list(srv120, med),funlabel=c("120-mo
> >> Survival","Median Survival"),fun.at=list(at.surv, at.med))
> >>
> >> I get a the following warning:
> >> Warning message:
> >> In approx(fu[s], xseq[s], fat) : collapsing to unique 'x' values
> >>
> >> However, a great nomogram is constructed.
> >>
> >>
> >> But then I try to stratify...
> >
> >>
> >> Stratified Nomogram:
> >>
> >> f<-cph(S~A+B+C+D+E+F+strat(H),x=T,y=T,surv=T,time.inc=10*12,method="b
> >> reslow")
> >>
> >> srv=Survival(f)
> >> surv.p <- function(lp) srv(10*12, lp, stratum="Hist=P") surv.f <-
> >> function(lp) srv(10*12, lp, stratum="Hist=F") surv.o <- function(lp)
> >> srv(10*12, lp, stratum="Hist=O")
> >> quant=Quantile(f)
> >> med.p <- function(lp) quant(.5, lp, stratum="Hist=P") med.f <-
> >> function(lp) quant(.5, lp, stratum="Hist=F") med.o <- function(lp)
> >> quant(.5, lp, stratum="Hist=O") nomogram(f, fun=list(surv.p, surv.f,
> >> surv.o, med.p, med.f, med.o),
> >> + funlabel=c("S(120|P)","S(120|F)","S(120|O)",
> >> + "med(P)","med(F)","med(O)"),
> >> + fun.at=list(c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9),
> >> + c(120,80,60,40,30,20,15,10,8,6,4,2,0)))
> >>
> >>
> >> the final nomogram only gives me a survival probability line for one
> >> of the
> >> 3 Hist categories "S(120|P)". It does show the letters "S(120|F)" but
> >> there is no survival probability line; there is nothing for the last
> >> category O, and no median risk at all.
> >
> > Those outputs seem consistent with the fact that stratification is not
> > computing separate models, but rather a pooled model. See Section
> > 19.1.7 of RMS.
>
> But you can think of stratification as using a different transformation
> for each stratum, and as long as you create a separate function for each
> level of the stratification variable, as Rene did, all should be well.
>
> >
> >> I considered the idea that I was exceeding some sort of space
> >> limitation, and tried to set total.sep.page=T, but it didn't change
> >> the output.
> >
> > Does a "median risk' exist when you stratify? You are allowing 3
> > separate survival functions to be created so that you estimate the
> > remaining parameters.
> > It's possible that you can extract information about them, but you may
> > be on your own about how to recombine them.
>
> Yes it exists, using the separate function approach.
>
> Rene if you can duplicate the problem with a simple simulated or real
> dataset and send that to me I can try to go through this step by step.
> It's probably a scaling, units of measurement, or extrapolation problem
> where the median is not defined. You can evaluate the created functions
> yourself a several settings to see if the results are reasonable and to
> learn where extrapolation is not possible because of truncated follow-up.
>
> Frank
>
> >
> >
> >> I get the following error message:
> >>
> >> Error in axis(sides[jj], at = scaled[jj], label = fat[jj], pos = y,
> >> cex.axis = cex.axis, :
> >> no locations are finite
> >>
> >>
> >> I would very much appreciate any assistance in this matter. Thank you
> >> very
> >> much.
> >>
> >> ~Renee Park
> >
> > David Winsemius, MD
> > Heritage Laboratories
> > West Hartford, CT
> >
> > ______________________________________________
> > 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-projectorg/posting-guide.html
> <http://www.R-project.org/posting-guide.html>
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
>
> --
> Frank E Harrell Jr Professor and Chair School of Medicine
> Department of Biostatistics Vanderbilt University
>
More information about the R-help
mailing list