# 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.
Also to speed up your R code note you can just call etimes = rexp(n, relhazs).
Note you do recover the actual estimates from using survreg(Surv(time, status) ~ x1 + x2, data=dat, dist='exp').
• Many thanks. One question: Do you know why survreg's output flips the sign of the coefficients? Call: survreg(formula = Surv(time, status) ~ x1 + x2, data = dat, dist = "exp") Coefficients: (Intercept) x1 x2 -0.002551373 -1.004660076 -2.993056428 Scale fixed at 1 Loglik(model)= -100274.5 Loglik(intercept only)= -655669.9 Chisq= 1110791 on 2 degrees of freedom, p= <2e-16 n= 100000 – Lei Huang Mar 18 at 20:47
• One more question/comment: What you said of the rounding issue sounds plausible; hazard ratios in this case do span many orders of magnitudes. But still, being between 10^-5 and 10^5, they are far from ranges where numerical underflow or overflow happens. In fact, when I tried two survival packages in Python (lifelines and scikit-survival), they are fine. Please see the additions to the question. Any idea why R's survival seems more brittle in handling this? – Lei Huang Mar 18 at 20:53
• You are correct about the rounding: coxph is detecting ties that don't exist (I assume if times are within some epsilon, it considers them a tie). You can turn this off with the following and recover the original coefficients: fit <- coxph(Surv(time, status) ~ x1+x2, data=dat, control = coxph.control(timefix = FALSE)) – Cliff AB Mar 18 at 20:58
• Thanks @CliffAB. Now it makes much more sense. Although I feel like the default behavior of survival's coxph should be timefix = FALSE, to avoid similar confusions in the future ;) – Lei Huang Mar 18 at 21:13
• For people's reference, the official documentation of coxph has mentioned this: A second issue is that of artificial ties due to floating-point imprecision. See the vignette on this topic for a full explanation or the timefix option in coxph.control. Users may need to add timefix=FALSE for simulated data sets. (rdrr.io/cran/survival/man/coxph.html#heading-9) – Lei Huang Mar 18 at 21:33