Under what conditions are HRs and RR consistent estimators of the same quantity?

Question

Assuming we have data for two groups from a time-to-event analysis with no censoring, there are two modeling approaches to comparing the group association with the risk of the event. Option one is a Cox proportional hazards model where a linear model is estimated having:

$$ \log \lambda (X|Z)= \log \lambda(X|Z=0) + \beta X$$

that is estimated by maximizing the partial likelihood

$$ \mathcal{L}(\beta; X, t) = \prod_{j=1}^K \frac{\lambda(X_j|Z_j)}{\sum_{\mathcal{l} \in\mathcal{R}_j\lambda(X_j|Z_\mathcal{l}) }} $$

More detail found here.

Alternately, one can consider the Poisson likelihood for an event and a log-linear model for a constant-in-time hazard function and an offset to designate the observation time, $t$ and each failure having the Poisson count variable equal to 1.

$$\log \lambda (X) = \alpha + \beta X - \log(t)$$

and $$ \mathcal{L}(\alpha, \beta, X) = \prod_{i=1}^n \lambda(X_i) \exp(-\lambda(X_i))$$

The question I have is under what conditions of the hazard function, or sample, is the $\beta$ that maximizes the Cox partial likelihood equal (or asymptotically equivalent to) to the $\beta$ that maximizes the Poisson likelihood? How do we prove this?

Answer 1

I'm not sure that there's a simple answer outside of the trivial case of a constant hazard.

As $\lambda$ is often used for the hazard function in a survival model, let's use $\mu_i = t_i\exp(\alpha)\exp(\beta X_i)$ to represent the Poisson mean for case $i$. The log-likelihood for the Poisson model (with all outcome values identically 1) is then:

$$\sum_i \log \mu_i - \mu_i=\sum_i \alpha + \beta X_i - t_i\exp(\alpha)\exp(\beta X_i)+ \log t_i.$$

The $\log t_i$ terms become 0 when differentiating with respect to the parameters to get the score function.

McCullagh and Nelder in Chapter 13 write the log-likelihood for a proportional hazards (PH) model with a specified baseline hazard in a way that parallels a Poisson model. With baseline hazard $h_0(t)$ and cumulative baseline hazard $H_0(t)$, then for case $i$ the cumulative hazard at its event time $t_i$ under PH is $\mu_i=H_0(t_i)\exp(\beta X_i)$. The log-likelihood without censored cases is:

$$ \sum_i \log \mu_i - \mu_i + \log \frac{h_0(t_i)}{H_0(t_i)}= \sum_i \log h_0(t_i) + \beta X_i - H_0(t_i)\exp(\beta X_i) . $$

For $h_0(t) = \exp(\alpha)$, this gives the same score function as the Poisson log-likelihood, as expected for an exponential survival model. The Cox partial likelihood approach can match the Poisson full-likelihood results for exponentially distributed data (Weibull with shape of 1) reasonably well.

> set.seed(101)
> expDF <- data.frame(time=c(rweibull(100,shape=1,scale=1),rweibull(100,shape=1,scale=2)),status=1,x=c(rep(0,100),rep(1,100)))
> expPois <- glm(status~x+offset(log(time)),data=expDF)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.6593     0.1406  11.801  < 2e-16 
x            -0.6000     0.1988  -3.018  0.00288  

> expCox <- coxph(Surv(time,status)~x,data=expDF)

     coef exp(coef) se(coef)      z Pr(>|z|)    
x -0.6416    0.5264   0.1494 -4.295 1.74e-05 

The ability to fit PH models with Poisson models having event-specific intercept terms has has been recognized for over 40 years: J Whitehead, "Fitting Cox's Regression Model to Survival Data using GLIM," Journal of the Royal Statistical Society. Series C (Applied Statistics) 29: 268-275 (1980).

Once you get beyond the exponential, however, you get into trouble. With a Weibull shape of 2 for simulated data (thus maintaining a PH interpretation amenable to Cox regression) instead of the exponential shape of 1:

> set.seed(101)
> weibDF <- data.frame(time=c(rweibull(100,shape=2,scale=1),rweibull(100,shape=2,scale=2)),status=1,x=c(rep(0,100),rep(1,100)))
> weibPois <- glm(status~x+offset(log(time)),data=weibDF)
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.32965    0.07030  18.914  < 2e-16
x           -0.64658    0.09942  -6.504 6.26e-10

> weibCox <- coxph(Surv(time,status)~x,data=weibDF)

     coef exp(coef) se(coef)      z Pr(>|z|)    
x -1.2927    0.2745   0.1680 -7.696  1.4e-14

The issue is that the Cox PH analysis doesn't take actual time into account, only the ordering of events. Any stretch or compression of the time scale between events will give the same coefficient estimates, regardless of how the underlying baseline hazard over actual time is thus changed. In contrast, the Poisson coefficient estimates depend explicitly on event times. This failure of agreement in a fairly simple case with non-constant hazard suggests that there might not be much to pursue in general, even though I can't rule out some other specific combinations of hazards and data that could provide agreement absent constant hazards in time.

The Poisson estimate agrees much better with the (negative of) the accelerated failure time (AFT) coefficient of interest (Weibull is the default for survreg()):

> weibSR <- survreg(Surv(time,status)~x,data=weibDF)
              Value Std. Error      z      p
(Intercept) -0.0370     0.0556  -0.67   0.51
x            0.6914     0.0766   9.03 <2e-16
Log(scale)  -0.6136     0.0557 -11.02 <2e-16

which I guess isn't that surprising, as the Poisson and AFT models (with all cases having events) substantially differ only in terms of the assumed error distributions of $\log t_i$ about the predicted values.