[R] Coxph not converging with continuous variable

[Keil, Alex]

The coxph function in R is not working for me when I use a continuous predictor in the model. Specifically, it fails to converge, even when bumping up the number of max iterations or setting reasonable initial values. The estimated Hazard ratio from the model is incorrect (verified by an AFT model). I've isolated it to the "x1" variable in the example below, which is log-normally distributed. The x1 here has extreme values, but I've been able to reproduce the problem intermittently with less extreme values. It seemed odd that I couldn't find this question asked anywhere, so I'm wondering if I'm just not seeing a mistake I've made.

Any help is appreciated to get this model to converge. Code, output, sessionInfo follows signature.

Alex Keil
UNC Chapel Hill


Here's an example:

library(survival)
set.seed(8182828)
#data
	N = 100000
	shape = 0.75
	hr = 3.5
	betax <- -log(hr)/(shape)
	ctime = 5

	z1 <- rbinom(N, 1, 0.5)
	z2 <- rbinom(N, 1, 0.5)
	x1 <- exp(rnorm(N, 0 + 2*z1, 1))
	wscale1 <- exp(0 + betax*x1 + z1 + z2 )
	time1 <- rweibull(N, shape, wscale1)
	t1 <- pmin(time1, rep(ctime, N))
	d1 <- 1*(t1 == time1)
	dat <- as.data.frame(cbind(t1, time1, d1, z1, z2, x1))
	 summary(dat)

#output
> summary(dat)
       t1               time1               d1               z1            z2               x1          
 Min.   :0.000000   Min.   : 0.0000   Min.   :0.0000   Min.   :0.0   Min.   :0.0000   Min.   :  0.0147  
 1st Qu.:0.000004   1st Qu.: 0.0000   1st Qu.:1.0000   1st Qu.:0.0   1st Qu.:0.0000   1st Qu.:  0.9552  
 Median :0.008400   Median : 0.0084   Median :1.0000   Median :0.5   Median :0.0000   Median :  2.7052  
 Mean   :0.295397   Mean   : 0.3260   Mean   :0.9906   Mean   :0.5   Mean   :0.4997   Mean   :  6.8948  
 3rd Qu.:0.178325   3rd Qu.: 0.1783   3rd Qu.:1.0000   3rd Qu.:1.0   3rd Qu.:1.0000   3rd Qu.:  7.7409  
 Max.   :5.000000   Max.   :52.8990   Max.   :1.0000   Max.   :1.0   Max.   :1.0000   Max.   :431.2779 

> exp((cox1 <- coxph(Surv(t1, d1)~  x1 + z1+ z2, ties="breslow"))[[1]]) #hrs
       x1        z1        z2 
3.3782387 0.4925040 0.4850214 
Warning message:
In fitter(X, Y, strats, offset, init, control, weights = weights,  :
  Ran out of iterations and did not converge


#accelerated failure time model confirms estimate is off from coxph

> aft1 <- survreg(Surv(t1, d1)~ x1 + z1 + z2, dist="weibull"); coefs <- aft1$coefficients; scale <- aft1$scale
> 	data.frame(psi = -coefs[2], scale, hr=exp(-coefs[2]*(1/scale))) #hr, scale is "weibull shape" in sas
        psi    scale       hr
x1 1.670406 1.333391 3.499958
#end output

> sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] splines   stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] survival_2.36-14

loaded via a namespace (and not attached):
[1] timereg_1.6-5