On Nov 20, 2011, at 6:34 PM, Paul Johnston wrote:
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).
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?
2. Should I calculate HR and CIs that "agree" with the LR or logrank P? How do I do that?
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. -- David.
Thank you, Paul
David Winsemius, MD West Hartford, CT ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
