[R] plotting spline term when frailty term is included

[David Winsemius]

> 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