muhaz()不返回基线危险率，而是返回包括 协变量对最终危险的贡献的危险函数。您可以像我对这段代码所做的那样划分这两个贡献。首先开始模拟具有固定危险率的数据集： $h_0$

library(survival)
set.seed(6)
n    <- 200
age  <- 50 + 12*rnorm(n)
sex  <- factor(sample(c('Male','Female'), n, rep=TRUE, prob=c(.6, .4)))
cens <- 15*runif(n)
h0   <- .02
h    <- h0*exp(.05*age-.12*(sex=='Female'))  # hazard
dt   <- -log(runif(n))/h                     # calculated survival time
e    <- ifelse(dt <= cens,1,0)               # if time is lower than censoring
                                             # time event = 1 else event = 0
dt   <- pmin(dt, cens)                       # the effective observation time
                                             # is the minimum between dt and cens
S    <- Surv(dt, e)                          # create the Survival function

f    <- coxph(S ~ age + sex)                 # define Cox model
cf   <- f$coefficients                       # value of coefficients
pval <- summary(f)$coefficients[,5]          # value of p-values

    # create the data.frame:
dataset     <- data.frame(cbind(age, sex, "time" = dt, "status" = e))   
dataset$sex <- factor(x = dataset$sex, labels = c("Male", "Female"))

之后，您可以使用 计算风险函数muhaz。我更喜欢通过设置 arguments 的值来创建一个更详细（即使计算密集型）的函数n.min.grid = 1000, n.est.grid = 2000：

fmuhaz              <- muhaz::muhaz(times = dataset$time, 
                                    delta = dataset$status, 
                                    n.min.grid = 1000, n.est.grid = 2000)
hfun                <- approxfun(x = fmuhaz$est.grid, y = fmuhaz$haz.est)  
dataset.order       <- dataset[with(dataset, order(time)),]      
dataset.order$BetaX <- exp( cf[1]*dataset.order$age + 
                            cf[2]*as.numeric(dataset.order$sex=="Male") ) 
dataset.order$h0    <- NA
for (l in 1:(nrow(dataset) - 9)) {
  dataset.order$h0[l] <- hfun(dataset.order$time[l]) /
                              mean(dataset.order$BetaX[l:nrow(dataset.order)])
}

plot(fmuhaz, lwd=2)          
  abline(h=h0, col="red", lty=2, lwd=2)
  lines(x=dataset.order$time[1:(nrow(dataset.order) - 9)], 
    y=dataset.order$h0[  1:(nrow(dataset.order) - 9)], lwd=2, col="blue")

根据时间协变量对数据集进行排序后，我创建了一个函数，该函数使用approxfun. 然后我计算有序数据集中每次的危险率值，并将该值除以其中是事件的数量，是事件的总数事件，换句话说，我只计算在计算时仍有风险的所有案例。当至少有 10 个以上的案例处于危险之中时，总和就会停止，因此有机会不让危险函数的值跳高。10的数字取自muhaz的帮助$\sum_{i}^{T}\beta'X$$i$$T$作者根据 和 给出了危险函数计算的fmuhaz$est.grid边界fmuhaz$haz.est。在您在输出中看到的最终绘图中，可以看到muhaz计算出的风险函数（黑线）、用于模拟数据集的原始风险函数（红色虚线）以及通过我解释的方法计算的基线风险函数（蓝线）。您可能会看到实际上蓝线和红线彼此靠近。