R help - Mar 2016 - plotting spline term when frailty term is included

Dear colleagues,

I'd very much appreciate your help in resolving a problem that I'm
having with plotting a spline term.

I have a Cox PH model including a smoothing spline and a frailty term as
follows:

   fit<-coxph(Surv(start,end,exit) ~ x + pspline(z) + frailty(a))

When I run a model without a frailty term, I use the following in order to plot
the splines:

   termplot(fit, term = 2, se = TRUE, ylab = "Log Hazard", rug=TRUE,
xlab = "z_name")

However, when the frailty term is included, it gives this error:

   Error in pred[, first] : subscript out of bounds

What am I doing wrong here? Or is it the case that termplot does not work when
splines and frailty are included?

A similar question was asked a number of years ago - as far as I can see it
didn't receive an answer (online, at least).
http://r.789695.n4.nabble.com/a-question-in-coxph-td908035.html

Best,
Conor Little



Conor Little
Postdoctoral Researcher
Department of Political Science, University of Copenhagen


	[[alternative HTML version deleted]]

> On Mar 1, 2016, at 2:59 AM, Conor Edward Little <cli at ifs.ku.dk>
wrote:
> 
> Dear colleagues,
> 
> I'd very much appreciate your help in resolving a problem that I'm
having with plotting a spline term.
> 
> I have a Cox PH model including a smoothing spline and a frailty term as
follows:
> 
>   fit<-coxph(Surv(start,end,exit) ~ x + pspline(z) + frailty(a))
> 
> When I run a model without a frailty term, I use the following in order to
plot the splines:
> 
>   termplot(fit, term = 2, se = TRUE, ylab = "Log Hazard",
rug=TRUE, xlab = "z_name")
> 
> However, when the frailty term is included, it gives this error:
> 
>   Error in pred[, first] : subscript out of bounds
I think this is due to the inability of `predict` to cope with the frailty model
when type="terms". Using the example on the ?frailty page as teh basis
for a reproducible example:
> predict(coxph(Surv(time, status) ~ age + frailty(inst, df=4), lung),
type="terms")
Error in predict.coxph.penal(coxph(Surv(time, status) ~ age + frailty(inst,  : 
  length of 'dimnames' [2] not equal to array extent
> 
> What am I doing wrong here? Or is it the case that termplot does not work
when splines and frailty are included?
> 
> A similar question was asked a number of years ago - as far as I can see it
didn't receive an answer (online, at least).
> http://r.789695.n4.nabble.com/a-question-in-coxph-td908035.html

The `?termplot` help-page says that there must be a `predict` method for the
object that accepts a "terms"- argument. Therneau advises the use of
the 'coxme'-package for mixed effects survival modeling and the
`predict.coxme` function does not have a "terms" argument.

Prediction with frailty: I think that the issue of prediction with frailty
models may be more complex than you understand. My level of understanding is
inadequate as well, so I tried searching for earlier answers. If Thomas Lumley
says it is "hard" then I take him at my word. Appears you would need
to do calculation separately for each value of the frailty "term" if
you stay with the coxph/frailty environment. Therneau has written: 
"Residuals methods for coxme would be an
important addition and is on my to-do list. (But as my wife would pointout, so
is a bathroom remodel and she isn't holding her breath.)".  So it
appears that the problem is more than just a simple two line bit of code.

Lumley's approach from 10+ years ago: 

http://markmail.org/search/?q=list%3Aorg.r-project.r-help+coxph+frailty+prediction+#query:list%3Aorg.r-project.r-help%20coxph%20frailty%20prediction%20+page:1+mid:qy7nxom3uknj2oop+state:results

http://markmail.org/search/?q=list%3Aorg.r-project.r-help+coxph+frailty+prediction+#query:list%3Aorg.r-project.r-help%20coxph%20frailty%20prediction%20+page:1+mid:nbcdjwfs3jhqjsji+state:results

> 
> Best,
> Conor Little
> 
> 
> 
> Conor Little
> Postdoctoral Researcher
> Department of Political Science, University of Copenhagen
> 
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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