how to fit Cox PH model in r (more than one factors)?

Question

I have a dataset for studying breast cancer patients. I wanted to fit a Cox proportional model (without considering the interaction term).

The variables contain age (<40 =1, 40~60=2,>60=3), predominant site (not middle=1, middle=2,unknown=9), maximum diameter (<2.5=1, 2.5~5.5=2), menopausal status (<2 year=1, >2 years=2,unknown=9), estrogen level ( neg=0, pos=1,unknown=9), progesterone levels (neg=0, pos=1,unknown=9) and w.censored(0=censored,1=not censored).

    'data.frame':   572 obs. of  6 variables:
 $ age       : Factor w/ 3 levels "1","2","3": 1 1 2 1 2 1 1 2 2 1 ...
 $ mepl.sts  : Factor w/ 3 levels "1","2","9": 1 3 2 1 2 1 1 2 1 1 ...
 $ pre.site  : Factor w/ 3 levels "1","2","9": 1 3 1 3 3 1 3 3 1 3 ...
 $ max.dia   : Factor w/ 4 levels "1","2","3","9": 1 2 4 3 3 3 3 3 2 2 ...
 $ es.level  : Factor w/ 3 levels "0","1","9": 3 3 3 2 2 1 1 1 1 1 ...
 $ prog.level: Factor w/ 3 levels "0","1","9": 3 3 3 1 1 1 1 1 1 1 ...

First I turned all these categorical variables into factors using as. factor. Then I did a Cox fit of all the variables in R and got the following results.

fit <- coxph(Surv(surv.day,w.cens) ~age + mepl.sts + pre.site+max.dia
               + es.level+prog.level ,data=bcnew)
summary(fit)
> summary(fit)
Call:
coxph(formula = Surv(surv.day, w.cens) ~ age + mepl.sts + pre.site + 
    max.dia + es.level + prog.level, data = bcnew)

  n= 572, number of events= 74 

               coef exp(coef) se(coef)      z Pr(>|z|)    
age2        -0.6059    0.5456   0.4408 -1.374 0.169287    
age3        -0.1771    0.8377   0.5457 -0.325 0.745463    
mepl.sts2    0.1884    1.2073   0.4145  0.455 0.649412    
mepl.sts9   -0.1904    0.8266   0.6041 -0.315 0.752580    
pre.site2    0.9555    2.5999   0.3594  2.659 0.007846 ** 
pre.site9    0.6220    1.8627   0.3260  1.908 0.056347 .  
max.dia2     1.0824    2.9518   0.3722  2.908 0.003632 ** 
max.dia3     1.9059    6.7256   0.4570  4.170 3.04e-05 ***
max.dia9    -0.8610    0.4227   0.6148 -1.400 0.161380    
es.level1   -0.6653    0.5141   0.3525 -1.887 0.059152 .  
es.level9   -1.1442    0.3185   0.3101 -3.690 0.000225 ***
prog.level1 -1.1256    0.3245   0.3921 -2.871 0.004095 ** 
prog.level9      NA        NA   0.0000     NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 

Concordance= 0.826  (se = 0.021 )
Likelihood ratio test= 103.8  on 12 df,   p=<2e-16
Wald test            = 94.27  on 12 df,   p=7e-15
Score (logrank) test = 140.6  on 12 df,   p=<2e-16

Since the p-value of age and menopausal status was > 0.1, the model was

h(t|x)=h0(t)exp(0.9555pre.site(2)+0.622pre.site(9)+1.0824max.dia(2)+1.9059max.dia(3) -0.6653es.level(1)-1.1442es.level(9)-1.1256*prog.level(1))

I don't know if this model is correct, but I think there is something strange about the results. It is common sense that a person's age would have an effect on survival time, but each age grouping gets a p-value much greater than 0.1.

By the way, if I have 20 variables grouped in my dataset, can I use the step(fit) procedure to get the final model?

Thank you very much!

Answer 1

One source of your problem is trying to fit too many coefficients on too small a data set. Although 572 sounds like a lot of cases, the information in a survival model is essentially determined by the number of events. You only have 74 events.

The usual rule of thumb to avoid overfitting in survival analysis is to have 10-20 events per coefficient that you are estimating, unless you are using some form of penalization. With 74 events, you should only be trying to fit 4 to 7 coefficients. You are fitting 12. That runs a couple of types of risk.

One is missing true associations with outcome, with high standard errors of coefficients due to small event numbers. That might be what's going on with age--adding more predictors to the model can diminish the apparent significance of other predictors that are associated with outcome.

The other is finding false associations with outcome that might happen to occur in this data set but wouldn't replicate in another--you might just fit noise in these data. In this case with so many predictors, I would guess that your reasonably high concordance of 0.82 with this model wouldn't be found on a new data set.

I'd recommend using Frank Harrell's course notes and book on Regression Modeling Strategies as a guide. You'll find there several other ways to improve your modeling, including:

• Don't set up separate categories for "unknown." Use imputation to estimate values in a way that allows the modeling to take that extra uncertainty into account. That's a well respected procedure that will help prevent bias.
• Don't bin continuous predictors like age into groups. Model them as continuous, and as flexibly as possible. For example, young people often have a more aggressive form of cancer (e.g., due to a genetic problem) than older individuals, so they die sooner after diagnosis. A U-shaped association of age with survival after diagnosis, as your data indicate, is quite possible. Your model should be able to handle that.
• Don't fit a model and then arbitrarily throw out predictors whose coefficients have p-values > 0.1. Don't focus so much on p-values at all. Avoid automated model selection, as you propose with step(fit).
• Use bootstrap resampling to validate and calibrate your model, checking how much overfitting might be involved.

Survival analysis is tricky. If you omit any predictor associated with outcome you can bias results for included predictors. But if you include more predictors than your number of events allow, you risk overfitting. The Harrell references should point you in the right direction.

Answer 2

If you know with fairly high certainty that a covariate has clinical relevance, you include it in the model even if it happens to fit badly on your training data. Essentially, your prior is strong enough that you would need large amounts of negative evidence to reject it.

Now why doesn't age seem significant in your model? I'm not familiar with the data so I can't tell. What I can tell you is that I would be much more methodical in developing the model.

A good first step is to study the individual covariates and their relationship to the outcome. Do that! Does age still look insignificant?