The question that arises in my mind, what benefit do you gain from the cox-model? You risk running into the proportional hazards problem and I am not sure to what your gain is. There is also the question if hazard ratio is what you want; the hazard ratio is a difficult concept while odds ratios are much more intuitive for most.
Alternative models
If you have one measurement per patient I would suggest doing a logistic regression with time as a non-linear variable. If you have multiple measurements per patient a random effects model is probably more useful. By having time in your model you the other variable's impact are independent of the time effect. Furthermore, you can easily test if some variables have an interesting interaction with the time, i.e. that some predictors speed up the healing process.
Your outcome
Furthermore, are you certain that you actually have healing as an outcome or is it just a threshold for a continuous score? Unless you are measuring infectious diseases or something similar, I would say that there are few things that are "healed", i.e. people have more or less symptoms after two months and it is a pity to throw away that information. I always try to explain this to my colleagues, the phrase: "you will tomorrow do surgery on an overweight patient" contains much less information than: "you will tomorrow do surgery on a patient with BMI 34". The latter is much more useful, contains more information and simplification doesn't always help.
Applying splines in R
I'm a big fan of splines and I've written several posts on the subject, here is one where you can see how to implement splines in R. You can also use piecewise linear. Bellow follows some code on how to apply these:
library(splines)
# The "true" disease process
# Observed days
days <- 1:70
# After 33 days everyone is healthy
# The score ranges from 0 to 10
disease <- ifelse(days < 33, days*.3, 10)
set.seed(1)
# Select a random sample from above
patients <- as.data.frame(cbind(days, disease)[sample(1:70, size=100, replace=TRUE),])
# Add noise to you'r observations
patients$observed_disease <- patients$disease + rnorm(100)
# Check the boundaries
patients$observed_disease[patients$observed_disease < 0] <- 0
patients$observed_disease[patients$observed_disease > 10] <- 10
fit <- lm(observed_disease ~ ns(days, 3), data=patients)
# The coefficients are uninterpreatable so you need to either
# 1. graph the results
# 2. contrast the results
pred_df <- as.data.frame(predict(fit, interval="confidence",
newdata=data.frame(days=1:70)))
library(ggplot2)
ggplot(patients, aes(x=days, y=observed_disease)) +
geom_point(color="#00000077") +
geom_ribbon(aes(x=days, ymax=upr, ymin=lwr, y=fit),
data=pred_df, fill="#00336644") +
geom_line(aes(x=days, y=fit), data=pred_df) +
theme_bw() +
theme(axis.title=element_text(size=24),
axis.text =element_text(size=14)) +
ylab("Disease score") +
xlab("Days")
The contrasting works as follows:
# The contrast function works best with
# the regression functions in the rms-package
library(rms)
ddist <- datadist(patients)
options(datadist="ddist")
fit.rms <- ols(observed_disease ~ rcs(days, 3), data=patients)
contrast(fit.rms,
a = list(days = c(10, 20, 30, 40, 50, 60)),
b = list(days = 30))
# Produces the following table:
# days Contrast S.E. Lower Upper t Pr(>|t|)
# 10 -5.506561 0.19106523 -5.8857731 -5.127350 -28.82 0
# 20 -2.472634 0.08149262 -2.6343748 -2.310894 -30.34 0
# * 30 0.000000 0.00000000 0.0000000 0.000000 NaN NaN
# * 40 1.333338 0.04217439 1.2496333 1.417042 31.61 0
# * 50 1.537336 0.09955526 1.3397462 1.734925 15.44 0
# * 60 1.138009 0.19349478 0.7539759 1.522043 5.88 0
#
# Redundant contrasts are denoted by *
#
# Error d.f.= 97
#
# Confidence intervals are 0.95 individual intervals
Piecewise linear regression may be of interest although it doesn't help you with your other variables, it is just a simplification:
library(segmented)
fit.segm <- lm(observed_disease ~ days, data=patients)
segm.out <- segmented(fit.segm, seg.Z=~days, psi=list(days=10))
summary(segm.out)
# Produces the following output:
#
# ***Regression Model with Segmented Relationship(s)***
#
# Call:
# segmented.lm(obj = fit.segm, seg.Z = ~days, psi = list(days = 10))
#
# Estimated Break-Point(s):
# Est. St.Err
# 32.5500 0.8944
#
# t value for the gap-variable(s) V: 0
#
# Meaningful coefficients of the linear terms:
# Estimat Std. Error t value Pr(>|t|)
# (Intercept) -0.15459 0.24450 -0.632 0.529
# days 0.30143 0.01178 25.594 <2e-16 ***
# U1.days -0.30204 0.01483 -20.364 NA
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.6923 on 96 degrees of freedom
# Multiple R-Squared: 0.9412, Adjusted R-squared: 0.9394
#
# Convergence attained in 5 iterations with relative change -4.622192e-16
# If you want to know the confidence interval for the breakpoint just do:
confint(segm.out)
# $days
# Est. CI(95%).l CI(95%).u
# 32.55 30.77 34.32
The above code should work with glm() if you are want to use logistic regression instead of linear regression.
