Skip to content

obtaining the values for the hazard function in a cox regression

2 messages · Bob Green, David Winsemius

Bob Green
# 
Hello ,

I am hoping for some advice regarding obtaining the values for the 
hazard function in a cox regression that I have undertaken. I have a 
model in the following form, analysed with the package survival (v. 
2.34-1) and a log-log plot obtained using Design (v. 2.1-2).

For two variables, the lines in the survival curves crossed. The 
statistician I been obtaining advice from (who does not use R) asked 
me to obtain the hazard function values. I am wanting to confirm 
whether basehaz is the correct command to obtain such values, to 
better understand what is occurring in the log plots.

Below is the code that I have been using.  Any assistance is appreciated,

regards

Bob



library(survival)
recidivismv <- read.csv("g://Chapsurv_v.csv",header=T)
cox.V <- coxph(Surv(intDaysUntilFVPO, Event_v) ~ intAgeAtMHCIndex + 
PRE + group + MHC + strGender, data = recidivismv)
summary (cox.V)


# produces a log-log plot

mvfit <- survfit (Surv(recidivismv$intDaysUntilFVPO, 
recidivismv$Event_v) ~ MHC, data = recidivismv)
neg.ll <- function(mvfit) -log(-log(mvfit))
library (Design)
survplot(mvfit, fun = neg.ll, conf = "none", logt=TRUE)

# hazard values
mhc2 <- coxph (Surv(recidivismv$intDaysUntilFVPO, 
recidivismv$Event_v) ~ MHC, data = recidivismv)
basehaz(mhc2)
David Winsemius
# 
It's not entirely clear from the help pages, but reading Thereau and  
Gramsch' "Modeling Survival Data" it appears to me that the basehaz  
function in the survival package returns the estimated baseline   
_cumulative_ hazard function H(t|mean covariates). Such a function  
will be monotonic upward. Search on "Nelson-Aalen estimator" for more  
detail on cumulative hazard functions. The literature on survival  
analysis seems ambiguous in that sometimes the cumulative hazard  
function is meant, and sometimes the instantaneous hazard function.

Clarify with your statistician what is needed. I found that a plot  
which reversed the default plotting  arrangement of the basehaz object  
was more helpful.

  with(basehaz.fit, plot(time, hazard, ylab = "Cumulative Hazard |  
Mean covariates") )