12
$\begingroup$

Can someone explain my Cox model to me in plain English?

I fitted the following Cox regression model to all of my data using the cph function. My data are saved in an object called Data. The variables w, x, and y are continuous; z is a factor of two levels. Time is measured in months. Some of my patients are missing data for variable z (NB: I have duly noted Dr. Harrell's suggestion, below, that I impute these values so as to avoid biasing my model, and will do so in the future).

> fit <- cph(formula = Surv(time, event) ~ w + x + y + z, data = Data, x = T, y = T, surv = T, time.inc = 12)

Cox Proportional Hazards Model
Frequencies of Missing Values Due to Each Variable
Surv(time, event)    w    x    y    z 
                0    0    0    0   14 

                Model Tests          Discrimination 
                                            Indexes        
Obs       152   LR chi2      8.33    R2       0.054    
Events     64   d.f.            4    g        0.437    
Center 0.7261   Pr(> chi2) 0.0803    gr       1.548    
                Score chi2   8.07                      
                Pr(> chi2) 0.0891                      

                   Coef    S.E.   Wald Z   Pr(>|Z|)
         w      -0.0133  0.0503    -0.26     0.7914  
         x      -0.0388  0.0351    -1.11     0.2679  
         y      -0.0363  0.0491    -0.74     0.4600  
         z=1     0.3208  0.2540     1.26     0.2067

I also tried to test the assumption of proportional hazards by using the cox.zph command, below, but do not know how to interpret its results. Putting plot() around the command gives an error message.

 cox.zph(fit, transform="km", global=TRUE)
            rho chisq      p
 w      -0.1125 1.312 0.2520
 x       0.0402 0.179 0.6725
 y       0.2349 4.527 0.0334
 z=1     0.0906 0.512 0.4742
 GLOBAL      NA 5.558 0.2347

First Problem

  • Can someone explain the results of the above output to me in plain English? I have a medical background and no formal training in statistics.

Second Problem

  • As suggested by Dr. Harrell, I would like to internally validate my model by performing 100 iterations of 10-fold cross-validation using the rms package (from what I understand, this would entail building 100 * 10 = 1000 different models and then asking them to predict the survival times of patients that they had never seen).

    I tried using the validate function, as shown.

    > v1 <- validate(fit, method="crossvalidation", B = 10, dxy=T)
    > v1
          index.orig training    test optimism index.corrected  n
    Dxy      -0.2542  -0.2578 -0.1356  -0.1223         -0.1320 10
    R2        0.0543   0.0565  0.1372  -0.0806          0.1350 10
    Slope     1.0000   1.0000  0.9107   0.0893          0.9107 10
    D         0.0122   0.0128  0.0404  -0.0276          0.0397 10
    U        -0.0033  -0.0038  0.0873  -0.0911          0.0878 10
    Q         0.0155   0.0166 -0.0470   0.0636         -0.0481 10
    g         0.4369   0.4424  0.6754  -0.2331          0.6700 10
    

    How do you perform the 100x resampling? I think my above code only performs the cross-validation once.

  • I then wanted to know how good my model was at prediction. I tried the following:

    > c_index <- abs(v1[1,5])/2 + 0.5
    > c_index
    [1] 0.565984
    

    Does this mean that my model is only very slightly better than flipping a coin?

Third Problem

Dr. Harrell points out that I have assumed linearity for the covariate effects, and that the number of events in my sample is just barely large enough to fit a reliable model if all covariate effects happen to be linear.

  • Does this mean that I should include some sort of interaction term in my model? If so, any advice as to what to put?
$\endgroup$
  • 2
    $\begingroup$ I made some major edits to my initial question, above, about three hours after I initially asked it. I have tried to follow Dr. Harrell's very helpful advice. I would still really appreciate it if someone would attempt to explain the above cph output to me in plain English, or point me to a reference that would do so. Dr. Harrell, thanks so much for your help so far! $\endgroup$ – Alexander Feb 17 '12 at 17:40
12
$\begingroup$

To get started, consider a few things. First, you are excluding too many observations with missing data and this will cause a bias. Consider multiple imputation. Second, there is a plot method for cox.zph which is useful in assessing proportional hazards. Third, you have assumed linearity for the covariate effects. Fourth, the number of events in your training sample is just barely large enough to fit a reliable model if all covariate effects happen to be linear (which is rare). And your test sample would have to have perhaps 400 events before it would yield a reliable assessment of prediction accuracy. It is not clear that you had enough data to split the data into two parts. Resampling validation (100 repeats of 10-fold cross-validation, or use the bootstrap) is a better solution. Both your original external validation (functions rcorr.cens and val.surv) and resampling internal validation (functions validate, calibrate) are implemented in the R rms package. Case studies for the rms package are found in my course notes at http://biostat.mc.vanderbilt.edu/rms (and I have a 3-day course on this in Nashville next month). Note that $2\times 2$ tables are not appropriate for use with continuous data.

$\endgroup$
  • 1
    $\begingroup$ Dr. Harrell, thanks a bunch for your comments. I tried typing plot(cox.zph(fit[[1]], transform="km", global=TRUE)), however, this yielded Error in plot.cox.zph(cox.zph(fit[[1]], transform = "km", global = TRUE)) : Spline fit is singular, try a smaller degrees of freedom. Am I calling this function incorrectly? $\endgroup$ – Alexander Feb 17 '12 at 14:45
  • 1
    $\begingroup$ If I combine my training and validation data, I have 166 observations with 75 events. As you suggested, I can impute the values of my variable "z" for my 13 observations who are missing that data, eliminating that problem. I like your suggestion of performing 100x 10-fold cross-validation. If you have time, I would really appreciate it if you could offer more of a concrete hint on how to use the rms packages to do this. In the mean time, I will continue to read through your website. I would love to take your course in the future. Unfortunately, I am in Europe at the moment, a bit too far away! $\endgroup$ – Alexander Feb 17 '12 at 15:14
  • 2
    $\begingroup$ For the singularity problem, do plot(cox.zph(...), df=2). See the case studies in the course notes for rms examples or install the package (which also needs the Hmisc package) and type these commands to bring up help files: ?cph ?validate.cph ?calibrate.cph $\endgroup$ – Frank Harrell Feb 17 '12 at 23:14
  • 2
    $\begingroup$ Thanks again for your help! This weekend I printed out your course notes and will spend this week reading them and working through all of the case studies. $\endgroup$ – Alexander Feb 20 '12 at 12:51
4
$\begingroup$

The output of the R cph function, based on a relevant example, is explained in this easy-to-follow paper by J. Fox.

I would strongly advise reading this paper if you have not already.

$\endgroup$
  • 1
    $\begingroup$ Can you describe how the paper suggests interpreting the cph output? $\endgroup$ – smillig Aug 5 '13 at 9:33
  • 2
    $\begingroup$ +1 Thanks for the reference and welcome to this site! It would be great if you could also give a succinct overview of the contents of the paper as we strive for answers that can stand on their own. $\endgroup$ – Gala Aug 5 '13 at 9:50
  • $\begingroup$ This link is not there anymore $\endgroup$ – Marcin Kosiński Jun 21 '15 at 13:21
  • 1
    $\begingroup$ The link socialsciences.mcmaster.ca/jfox/Books/Companion-1E/… is currently working and looks to be the same paper referenced in this answer. $\endgroup$ – dnlbrky Mar 26 at 13:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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