What does muhaz return? A pretty basic question.  I have read somewhere that the muhaz function in muhaz package will return the baseline hazard rate for COX model.  The muhaz document states that it "Estimates the hazard function from right-censored data using kernel-based methods"  So it looks to me that muhaz estimates the complete hazard function, not the COX model baseline hazard rate.  Can someone clarify?
Thanks  
 A: muhaz() doesn't return the baseline hazard rate, but the hazard function including  the contribution of covariates to the final hazard. You can divide the two contribution as I do with this code. First start simulating a dataset with a fixed hazard rate $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"))

Afterwards you can calculate the hazard function by using muhaz. I prefer to create a more detailed (even if computationally intensive) function by setting the value of 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")


After ordering the dataset according the time covariate I create the function that gives the values of the hazard rate by using approxfun. Then I calculate the values of the hazard rate for each time in the dataset ordered and divide the value by the $\sum_{i}^{T}\beta'X$ where $i$ is the number of the event and $T$ is the total number of the event, in other words I count only all cases remaining at risk at the moment of calculation. The sum stops when there are at least 10 more cases at risk, so giving the chance not to give high jumps in the value of hazard function. The number of 10 is taken from the help of muhaz where the authors give the boundaries of the calculation of the hazard function in terms of fmuhaz$est.grid and fmuhaz$haz.est. In the final plot you see in output it's possible to see the muhaz calculated hazard function (the black line), the original hazard function used for simulating the dataset (red dash line), and the baseline hazard function calculated by the method I explained (the blue line). You may see that actually the blue and red line lie close each other.
