I am trying to use the segmented package in R to identify breakpoints/change points in a Cox PH model. Specifically, I'm interested in finding the optimal number of breakpoints and then estimating the breakpoints.
Here's a working example of my code:
library(survival)
library(segmented)
cox_model <- coxph(Surv(time, status)~age+pat.karno+meal.cal, data=lung)
os<-segmented.default(cox_model, ~age, psi=c(60, 70), n.boot=50) # estimate the breakpoint in the age effect
By providing more values for psi (e.g. c(60, 65, 70), I notice that more breakpoints are identified. Is there a way to allow the function to automatically determine the optimal number of breakpoints?
Edit:
Apparently the selgmented function can help with this, however I am not sure how to apply it on the Cox model. Find my code below:
cox_model <- coxph(Surv(time, status)~age+pat.karno+meal.cal, data=lung)
os <-selgmented(cox_model, Kmax=10, type="bic") #BIC-based selection
It is giving me the following error:
Error in bic.values[1] <- BIC.f(olm) + 1 : replacement has length zero
In addition: Warning messages:
1: In min(Z) : no non-missing arguments to min; returning Inf
2: In max(Z) : no non-missing arguments to max; returning -Inf
Am I using it wrong?