Is it possible to specify baseline hazard in coxph?

Question

I am running simulations and I need to run a Cox model with a certain baseline hazard. The data are generated with a constant baseline hazard $A$ and I need to run coxph using this data but instead of estimating baseline hazard using Breslow formula, I need R to use baseline hazard=$B$). Is there a way to do it?

Answer 1

You're just doing exponential survival, e.g. parametric survival. See ?survreg in the survival package. If you specify dist='exponential' as an argument, it will assume a constant hazard function and estimate that baseline hazard along with the hazard ratios. It will not, for instance, take a user supplied baseline hazard and estimate the hazard ratios of the regression model.

As you know, when doing Cox regression, you make no assumption about the functional form of the baseline hazard, it is not estimated but factors out of the partial likelihood. So if the hazard is indeed constant, the Cox model will still estimate the right hazard ratio. What you lose by doing this is a bit of efficiency. In very small sample sizes, this is a legitimate concern.

If you want to estimate the HR with maximum likelihood, you can tweak the equations to give any baseline hazard estimate you want, even the "wrong" one ($B$ as you call it). Parametric survival models are just maximum likelihood. To do that, you need to modify the likelihood equations and come up with an estimating equation for the survival model. It shouldn't be too hard since you simply need to combine the exponential density functions dexp and pexp (for censored obs) along with nlm.

Example below:

library(survival)
set.seed(1)

## sample size
n <- 1000

## binary exposure (prevalence)
p <- 0.1 
x <- rbinom(n, 1, p)

## baseline hazard / hazard ratio
a <- 0.05
h <- 1.30

## generate survival times
t <- rexp(n, rate=exp(-log(a) + log(h)*x))

## check with cox model
coxph(Surv(t) ~ x)

## parametric exponential survival model (comes from "rate" based model so opposite)
survreg(Surv(t) ~ x, dist='exp')

## unconstrained ML, gives same results as above
negLogLik <- function(param) {
  -sum(dexp(t, exp(-param[1] + param[2]*x), log=TRUE))
}
nlm(negLogLik, c(0,0))

## now with the "right" BL hazard
negLogLikBLH <- function(param) {
  -sum(dexp(t, exp(-log(a) + param[2]*x), log=TRUE))
}
exp(nlm(negLogLikBLH, c(0,0))$estimate)

## now with the "wrong" BL hazard
b <- 0.10
negLogLikBLH2 <- function(param) {
  -sum(dexp(t, exp(-log(b) + param[2]*x), log=TRUE))
}
nlm(negLogLikBLH2, c(0,0))