1
$\begingroup$

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]
$\endgroup$
2
$\begingroup$

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').

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Lei Huang Mar 18 at 20:47
  • $\begingroup$ 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? $\endgroup$ – Lei Huang Mar 18 at 20:53
  • 2
    $\begingroup$ 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)) $\endgroup$ – Cliff AB Mar 18 at 20:58
  • $\begingroup$ 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 ;) $\endgroup$ – Lei Huang Mar 18 at 21:13
  • 1
    $\begingroup$ 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) $\endgroup$ – Lei Huang Mar 18 at 21:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.