I would really appreciate if you could help me with your expertise:

I have a study in which I evaluate patients for healing after discharge. So I have two group, healed and non-healed patients.

The question I want to ask is how to analyze independent predictors of healing (or non healing) within a defined time interval like within two months of the insult. The time for follow exams is variable but most follow ups were in the period of two months of the insult with some extending up to 6 months of insult.

I thought of using a cox model, where the censored data will be healed as non-healed patients which were examined after the period of interest were non-healed before that for sure but this is not the case for the healed ones. So wouldn't the Cox be reflective for my period of interest since most of the data lies within this period?

Thank you very much.

  • $\begingroup$ Yes I am going to use Cox regression for that most likely. But the risk factors identified will reflect which period of time? I think it will reflect the median time in the data (40 days)? Because I want to answer a question which is predctors of healing (or non-healing) within 2 months of discharge. $\endgroup$
    – SafSaf
    Commented Apr 3, 2014 at 13:13

1 Answer 1


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:


# 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)

# 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", 
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() + 
    axis.text =element_text(size=14)) + 
  ylab("Disease score") +

GGplot example

The contrasting works as follows:

# The contrast function works best with
# the regression functions in the rms-package
ddist <- datadist(patients)
fit.rms <- ols(observed_disease ~ rcs(days, 3), data=patients)
         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:

fit.segm <- lm(observed_disease ~ days, data=patients)
segm.out <- segmented(fit.segm, seg.Z=~days, psi=list(days=10))
# 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:
# $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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.