# R package survival: the estimated parameters are off

I simulated some survival data and used the R package survival's function coxph to estimate the parameters, but they are different from the parameter values I used for the simulation. What have I done wrong?

library(mvtnorm)
library(survival)

n <- 100000
set.seed(0)

# add some correlation, but the observation is similar for no correlation
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
etimes <- c()
for(i in 1:n){
# assume baseline hazard is a constant function at 1
# so the survival times are simply exp distributed
etimes[i] <- rexp(n=1, rate=relhazs[i])
}

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

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

fit <- coxph(Surv(time, status) ~ x1+x2, data=dat)
summary(fit)

Output:

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

n= 100000, number of events= 1e+05

coef exp(coef)  se(coef)     z Pr(>|z|)
x1  0.799999  2.225538  0.004016 199.2   <2e-16 ***
x2  2.391868 10.933898  0.006161 388.2   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

exp(coef) exp(-coef) lower .95 upper .95
x1     2.226    0.44933     2.208     2.243
x2    10.934    0.09146    10.803    11.067

Concordance= 0.932  (se = 0.001 )
Rsquare= 0.883   (max possible= 1 )
Likelihood ratio test= 214574  on 2 df,   p=<2e-16
Wald test            = 181957  on 2 df,   p=<2e-16
Score (logrank) test = 159179  on 2 df,   p=<2e-16

I expect to see the coef in the output to be close to 1 and 3 for x1 and x2, respectively. I did a similar simulation with Python's lifelines, and was able to recover the parameter values used in the simulation.

Thanks for any help!

Addition: Python codes that work using packages lifelines or scikit-survival:

import numpy as np
import pandas as pd
from lifelines import CoxPHFitter
from sksurv.linear_model import CoxPHSurvivalAnalysis

n = 10000
b1, b2 = 1, 3
np.random.seed(0)

X = pd.DataFrame(np.random.multivariate_normal([0,0], [[1,0.5],[0.5,1]], size=n),
columns=['x1','x2'])
x1, x2 = X.x1, X.x2
r = np.exp(b1*x1 + b2*x2)

etimes = [np.random.exponential(scale=1/r_) for r_ in r]

df = pd.DataFrame({'time':etimes, 'status':[1]*len(etimes), 'x1':x1, 'x2':x2})
y = np.array([(True, etime) for etime in etimes], dtype=[('status', '?'), ('time', '<f8')])

CoxPHFitter().fit(df, duration_col='time', event_col='status').print_summary()

print('\n', CoxPHSurvivalAnalysis().fit(X, y).coef_)

Output:

<lifelines.CoxPHFitter: fitted with 10000 observations, 0 censored>
duration col = 'time'
event col = 'status'
number of subjects = 10000
number of events = 10000
log-likelihood = -70242.56
time fit was run = 2019-03-18 20:38:54 UTC

---
coef  exp(coef)  se(coef)      z      p  -log2(p)  lower 0.95  upper 0.95
x1  1.02       2.78      0.01  72.01 <0.005       inf        0.99        1.05
x2  3.02      20.55      0.03 114.18 <0.005       inf        2.97        3.07
---
Concordance = 0.90
Log-likelihood ratio test = 23732.74 on 2 df, -log2(p)=inf

[1.02114613 3.02275361]

The issue is rounding. Log hazard ratios of 1 and 3 in normally distributed data ranging from -4 to 4 on average results in huge range of survival times, where half are less than 1, but the 99th percentile is 3682. Setting $$b_1$$ and $$b_2$$ to 1/10th their values reliably generates consistent estimates.