Bias in parameter estimates for Cox proportional hazard model when covariates are collinear

Question

For linear regression, if $y$ actually depends on two positively correlated covariates $x_1$ and $x_2$ (we can call it the true model), and if we only include one covariate, say $x_1$, in the regression model (we can call it the working model), its coefficient $\beta_1$ will be overestimated. This makes intuitive sense, because now $\beta_1$ represents the effect of both $x_1$ and $x_2$. One can in fact derive that $\tilde{\beta}_1 = \beta_1 + \rho \beta_2$, where $\tilde{\beta}_1$ is the apparent coefficient of $x_1$ in the working model, $\beta_1$ and $\beta_2$ are the actual coefficients of $x_1$ and $x_2$ in the true model, and $\rho$ is the correlation coefficient between $x_1$ and $x_2$. (The derivation is appended at the end.)

Now I am interested in the same question for Cox proportional hazard model. To my surprise, I observe that when the true model has positively-correlated $x_1$ and $x_2$, and the working model has only $x_1$, $\beta_1$ is in fact underestimated, at least when estimated using Cox's partial likelihood method. Here are my simulation codes with some explanatory comments.

library(mvtnorm)
library(survival)

n <- 100000
set.seed(0)

# set the correlation coefficient to be 0.5
sigma <- matrix(c(1,0.5,0.5,1), ncol=2)
X <- rmvnorm(n=n, mean=c(0,0), sigma=sigma)

x1 <- X[,1]
x2 <- X[,2]

b1 <- 1
b2 <- 3

# relative hazards
relhazs <- exp(b1*x1 + b2*x2)

# event times
# assume baseline hazard is a constant function at 1
# so the survival times are simply exp distributed
etimes <- rexp(n, relhazs)

# assume no censorship for simplicity
status <- rep(1, n)

dat <- data.frame(id=1:n,
                  time=etimes,
                  status=status,
                  x1=x1,
                  x2=x2)

Output:

Call:
coxph(formula = Surv(time, status) ~ x1 + x2, data = dat, control = coxph.control(timefix = FALSE))

        coef exp(coef)  se(coef)     z      p
x1  1.004289  2.729965  0.004424 227.0 <2e-16
x2  2.991976 19.925012  0.008271 361.7 <2e-16

Likelihood ratio test=236231  on 2 df, p=< 2.2e-16
n= 100000, number of events= 1e+05

Call:
coxph(formula = Surv(time, status) ~ x1, data = dat, control = coxph.control(timefix = FALSE))

       coef exp(coef) se(coef)     z      p
x1 0.837905  2.311519 0.003825 219.1 <2e-16

Likelihood ratio test=49198  on 1 df, p=< 2.2e-16
n= 100000, number of events= 1e+05

My questions:

1. Why the underestimation?
2. Can we possibly derive some analytical relationship between $\tilde{\beta}_1$, $\beta_1$, $\beta_2$ and $\rho$ in this case, just like what we did for linear regression, even for some simple baseline distribution such as $\text{Exp}(1)$ in my simulation?
3. If I estimate the parameters using a parametric model, then $\beta_1$ is indeed overestimated (please see below). Why the difference between Cox's semi-parametric partial-likelihood-based estimation and parametric full-likelihood-based estimation?

Call:
survreg(formula = Surv(time, status) ~ x1 + x2 + 0, data = dat, 
    dist = "exp")

Coefficients:
       x1        x2 
-1.004653 -2.993058 

Scale fixed at 1 

Loglik(model)= -100274.9   Loglik(intercept only)= -655669.9
    Chisq= 1110790 on 1 degrees of freedom, p= <2e-16 
n= 100000 

Call:
survreg(formula = Surv(time, status) ~ x1 + 0, data = dat, dist = "exp")

Coefficients:
      x1 
-2.50526 

Scale fixed at 1 

Loglik(model)= -1302576   Loglik(intercept only)= -655669.9
n= 100000 

---

Appendix:

For linear regression, suppose the true model is $y=\beta_1 x_1 + \beta_2 x_2 + \epsilon$, and the working model is $y=\tilde{\beta}_1 x_1 + \tilde{\epsilon}$. Without the loss of generality, assume $x_1$ and $x_2$ are standardized, so $x_2=\rho x_1$ where $\rho$ is the correlation coefficient. Equating the $y$ in the two equations and plugging in $x_2=\rho x_1$ give $\tilde{\beta}_1 = \beta_1 + \rho \beta_2$.

Answer 1

The Cox model effectively uses weighted least squares, with weights that are functions of the data and of $\hat{\beta}$. It also has the issue that the hazard ratio is not collapsible, meaning that omitting even a completely orthogonal predictor that is related to the outcome will cause a distortion in the remaining $\beta$s.

Answer 2

The situation is worse, in fact, and I think more comical. Try setting the correlation coefficient to 0, and run the same simulation. Surely, we should have less bias, maybe even close to no bias, right? No! We actually observe even more bias (~0.84 vs ~0.28). @Frank Harrell mentioned the non-collapsibility of the hazard ratio, but there is another bias occurring, a causal bias (this may be the same thing as the non-collapsibility, but an alternative perspective).

This causal bias occurs because the likelihood for the Cox model computes survival at time $t+1$ given survival at time $t$. This is fine if we can block all back door paths, but by omitting any risk factor (not even confounders), we fail to block the back door paths. Below I explain why. See the causal DAG below:

More specifically, when we wish to compute the effect of X1 on S_{t+1}, we (implicitly) condition on S_t. In the DAG, S_t is a collider, and this should worry us. For if there is another variable, like X2 that also effects survival, then conditioning on S_t opens up the back door path S_{t+1} <- X2 -> S_t <- X1. So, generally, for any t > 0, we have the effect of X1 being biased in an unpredictable way by any and all other risk factors,

A terrific paper on this, and where I learned this fact, is here: https://core.ac.uk/download/pdf/144149054.pdf