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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
