I am running a Cox proportional hazards model in R and need to include a time-by-covariate interaction (i.e., f(time) X covariate) in the model due to the violation of the proportional hazards assumption. My ultimate goal is to create a graph that shows how the hazard ratio of the covariate (which is binary) changes over time, using the coefficient on the main effect of covariate and the interaction term on time and covariate. I would like to use the restricted cubic splines for the time variable since time is non-linear in my case. I have been searching online how to do this in R but there are largely two things I'm still struggling with.
First, I use this code to run the Cox regression including the covariate-by-time interaction when time is a non-linear function estimated using the restricted cubic splines. However, the output shows coefficients on the interaction between the reference group of the 'key_predictor' (the group that equals to 0) and 'functime', rather than the interaction between the group that equals to 1 in the key binary predictor and 'functime'. This is solved when I include the main effect of 'functime' but I'm not sure it makes sense to include the main effects of time since the coefficients are NA's when I include their main effects. How can I get correct coefficients on the interaction terms without adding the main effect of time?
data <- data %>% mutate(time = tstop, functime = rcs(time, iknots=c(6,8,12), bknots=c(2,14))
fit <- coxph(Surv(tstart, tstop, endpt) ~ key_predictor + key_predictor:functime + other_controls, data= data)
Second, I used the following code to create a graph showing hazard ratio of 'key_predictor" over time (=months, ranging from 1 to 15). There is an error message saying "error in xy.coords(x, y, xlabel, ylabel, log) : 'x' and 'y' lengths differ" after I run those lines, but I don't know if the code below is a correct way to create such a graph in the first place. How can I create a graph showing the hazard ratio of a predictor (on the y-axis) over time (on the x-axis) taking into account the fact that time is modeled as restricted cubic splines?
t <- seq(1, 15, by = 1)
plot(x = t, y = coef(fit)["key_predictor"] + coef(fit)["key_predictor:functime"] * rcs(t, iknots=c(6,8,12), bknots=c(2,14)),
ylab = "Beta(t) for the key predictor", xlim = c(0, 16))