Skip to content

Cox proportional hazards confidence intervals

3 messages · Paul Johnston, David Winsemius, Peter Dalgaard

Paul Johnston
# 
I am calculating cox propotional hazards models with the coxph
function from the survival package.  My data relates to failure of
various types of endovascular interventions.  I can successfully
obtain the LR, Wald, and Score test p-values from the coxph.object, as
well as the hazard ratio as follows:

formula.obj = Surv(days, status) ~ type
coxph.model = coxph(formula.obj, df)
fit = summary(coxph.model)
hazard.ratio = fit$conf.int[1]
lower95 = fit$conf.int[3]
upper95 = fit$conf.int[4]
logrank.p.value = fit$sctest[3]
wald.p.value = fit$waldtest[3]
lr.p.value = fit$logtest[3]

I had intended to report logrank P values with the hazard ratio and CI
obtained from this function.  In one case the P was 0.04 yet the CI
crossed one, which confused me, and certainly will raise questions by
reviewers.  In retrospect I can see that the CI calculated by coxph is
intimately related to the Wald p-value (which in this specific case
was 0.06), so this would appear to me not a good strategy for
reporting my results (mixing the logrank test with the HR and CIs from
coxph).

I can report the Wald p-values instead, but I have read that the Wald
test is inferior to the score test or LR test.  My questions for
survival analysis jockeys out there / TT:

1. Should I just stop here and use the wald.p.value?  This appears to
be what Stata does with the stcox function (albeit Breslow method).

2. Should I calculate HR and CIs that "agree" with the LR or logrank
P?  How do I do that?

Thank you,
Paul
David Winsemius
#
On Nov 20, 2011, at 6:34 PM, Paul Johnston wrote:

            

      
      
      
    

        
I don't understand two things: Why would your report the inferior  
result, and I suppose  I also wonder why does it make that much  
difference? The estimate is what it is and a p-value of .04 is not  
that different from one of .06. Or are we dealing with religious  
beliefs here?

            

      
      
      
    

        
Therneau and Grambsch show how to calculate profile likelihood curves  
that can be used to generate confidence intervals on pages 57-59 of  
"Modeling Survival Data". This "survival analysis jockey" considers  
that book an essential reference. They basically use the offset  
capacity to construct 50 likelihoods around the estimate for one  
particular variable in a more complete model and then show where the  
97.5th and 0.025th percentile points are for an beta estimate based on  
a chi-square distribution for these log-likelihoods.

Further code not possible in the absence of the complete formula.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

      
On Nov 21, 2011, at 05:50 , David Winsemius wrote:

        

            

      
      
      
    

        
Well, it is an annoying inconsistency, but one that happens all over the place. 

Consider the single binomial sample: Most applied statistics textbooks teach the students to calculate the error margin for a CI as 1.96*sqrt(phat*(1-phat)/n), but the cutoff for testing p=p0 uses 1.96*sqrt(p0*(1-p0)/n). This will (and does) give rise to situations where the test and the CI disagrees. 

It is fixable by using the error margins at the upper and lower confidence limits, but it leads to nonlinear equations that you'd rather not inflict on students. (For the binomial sample it's a quadratic equation and prop.test actually uses it.)

I'd just leave it, and, if anyone complains, put in a note to the effect that "the slight inconsistency between p-value and CI is due to different large-sample approximations".