On 03/02/2016 05:00 AM, r-help-requ...@r-project.org wrote:
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?
There are 3 parts to the answer.
1. The first is a warning: wrt to mixed effects Cox models, I shifted my energy to coxme
over 10 years ago. The "penalized add on to coxph" approach of the frailty function was
an ok first pass, but is just too limited for any but the simplest models. I'm unlikely
fix issues, since there are others much higher on my priority list.
2. As Dave W. pointed out, the key issue is that predict(type='terms') does not work
with for a model with two penalized terms, when one is frailty and the other pspline.
Termplot depends on predict.
3. Again, as Dave W pointed out, the whole issue of what the "correct" answer would be
gets much more complicated when one adds random effects to the mix; some things have not
done because it is not clear where to go. (Survival curves for a mixed effects model only
recently got added to my "todo" list, even though it has been on the wish list forever,
because I finally have a notion of what a good approach would be.)
In your case I'd advise an end run: fit the model using ns() instead of pspline. I like
smoothing splines better than regression splines, but the fact is that for most data sets
they result in nearly identical answers.
Terry T
______________________________________________
R-help@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.
