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] （Score (logrank) test）
wald.p.value = fit$waldtest[3] （Wald test）
lr.p.value = fit$logtest[3] ---- Likelihood ratio test (优先）
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,