R help - May 2013 - Post-fit survival curve for shared frailty model

Hi,

I have fitted a shared frailty model using coxph in the survival
package with 'center' as random effect. Now I'm trying to obtain the
conditional survival curve for a new test case that comes from one of
the specified centers in the model.

The data set on which the model is fitted contains 1000 individuals
coming from 10 different centers, with 100 individuals in each center.

The fitted model has the following form:

model.frail = coxph(Surv(time,status) ~ x1 + frailty(center), data = d)

My new test case belongs to center no. 5 and has x1-value equal to 0.
My attempt to obtain the conditional survival function:

sfit = survfit(model.frail, newdata=data.frame(x1 = 0, center = 5))

This generates the error:
Error in survfit.coxph(model.frail, newdata = data.frame(x1 = 0,  :
  Newdata cannot be used when a model has sparse frailty terms

Am I specifying the newdata argument in the wrong way? If so, how
should it be specified?

I found an earlier note on the issue of frailty and survival curves by
Terry indicating that it is indeed possible to obtain the conditional
survival curve:

BEGINQUOTE

With respect to Cox models + frailty, and post-fit survival curves.

 1. There are two possible survival curves, the conditional curve
where we identify which center a subject comes from, and the marginal
curve where we have integrated out center and give survival for an
"unspecified" individual.  I find the first more useful.  More
importantly to your case, the survival package currently has no code
to calculate the second of these.

 2. When the number of centers is large the coxph code may have used a
sparse approximation to the variance matrix, for speed reasons.  In
this particular case one cannot use the "newdata" argument.  The
reason is entirely practical --- the code turned out to be very hard
to write. The need for this comes up very rarely, and the work around
is to use
   coxph(....... + frailty(center, sparse=1000, ....)
where we set the "sparse computation" threshold to be some number
larger than the number of centers, i.e., force non-sparse computation.

Terry Therneau

ENDQUOTE

Thanks in advance.

/Ulf

On May 31, 2013, at 6:01 AM, Ulf Sandvik wrote:
> Hi,
> 
> I have fitted a shared frailty model using coxph in the survival
> package with 'center' as random effect. Now I'm trying to
obtain the
> conditional survival curve for a new test case that comes from one of
> the specified centers in the model.
> 
> The data set on which the model is fitted contains 1000 individuals
> coming from 10 different centers, with 100 individuals in each center.
> 
> The fitted model has the following form:
> 
> model.frail = coxph(Surv(time,status) ~ x1 + frailty(center), data = d)
As I read the cited note below, this should have been:

 model.frail = coxph(Surv(time,status) ~ x1 + frailty(center, sparse=200), data
= d)

(Where 200 is arbitrary but larger than the 10 unique values. It's described
in ?frailty. Looking at the code I think it gets passed to
'frailty.gamma' in the default case.)

(There is a note on the ?frailty page: "The coxme package has superseded
this method. It is faster, more stable, and more flexible.")

-- 
David.
> 
> My new test case belongs to center no. 5 and has x1-value equal to 0.
> My attempt to obtain the conditional survival function:
> 
> sfit = survfit(model.frail, newdata=data.frame(x1 = 0, center = 5))
> 
> This generates the error:
> Error in survfit.coxph(model.frail, newdata = data.frame(x1 = 0,  :
>  Newdata cannot be used when a model has sparse frailty terms
> 
> Am I specifying the newdata argument in the wrong way? If so, how
> should it be specified?
> 
> I found an earlier note on the issue of frailty and survival curves by
> Terry indicating that it is indeed possible to obtain the conditional
> survival curve:
> 
> BEGINQUOTE
> 
> With respect to Cox models + frailty, and post-fit survival curves.
> 
> 1. There are two possible survival curves, the conditional curve
> where we identify which center a subject comes from, and the marginal
> curve where we have integrated out center and give survival for an
> "unspecified" individual.  I find the first more useful.  More
> importantly to your case, the survival package currently has no code
> to calculate the second of these.
> 
> 2. When the number of centers is large the coxph code may have used a
> sparse approximation to the variance matrix, for speed reasons.  In
> this particular case one cannot use the "newdata" argument.  The
> reason is entirely practical --- the code turned out to be very hard
> to write. The need for this comes up very rarely, and the work around
> is to use
>   coxph(....... + frailty(center, sparse=1000, ....)
> where we set the "sparse computation" threshold to be some number
> larger than the number of centers, i.e., force non-sparse computation.
> 
> Terry Therneau
> 
> ENDQUOTE
> 
> Thanks in advance.
> 
> /Ulf
> 
> ______________________________________________
> 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-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
David Winsemius
Alameda, CA, USA